Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Chem., 04 September 2025

Sec. Green and Sustainable Chemistry

Volume 13 - 2025 | https://doi.org/10.3389/fchem.2025.1656180

This article is part of the Research TopicAdvances in Chemical Looping Conversion of FuelsView all articles

Multiscale reaction model coupling dual-site microkinetics with bulk diffusion and CFD–DEM for a perovskite oxygen carrier

Ruiwen WangRuiwen Wang1Zhenshan Li
Zhenshan Li1*Lei LiuLei Liu2
  • 1Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Department of Energy and Power Engineering, Tsinghua University, Beijing, China
  • 2Hunan Engineering Research Center of Clean and Low-Carbon Energy Technology, School of Energy Science and Engineering, Central South University, Changsha, China

Non-catalytic heterogeneous reactions in fluidized beds involve physical and chemical processes spanning across the atom, surface, grain, particle, and reactor scales. However, a multiscale modeling framework covering all scales has not been fulfilled due to the incomplete coupling strategies. This study develops a multiscale model coupling all five scales. The elementary reaction path is derived from first-principles calculation, which is applied to a dual-site mean-field microkinetics describing the states of active site pairs; bulk-phase ion diffusion is treated by a lumped parameter method considering the asymmetrical effects of different site types. The intrinsic reaction kinetics is coupled with intraparticle gas diffusion and fluidization computed via CFD–DEM; experimental validation is conducted on a micro-fluidized-bed thermogravimetric analyzer measuring the solid conversion. The model is applied to the reduction of CaMn0.375Ti0.5Fe0.125O3−δ by H2 at designed gas concentrations and temperatures, revealing the effects of parameters from all scales on the overall reaction kinetics. The developed multiscale framework can be further adopted in other heterogeneous reactions with determined solid microstructures.

1 Introduction

Non-catalytic heterogeneous reactions between gases and solid crystals are involved in chemical looping (Ishida and Jin, 1994; Lyngfelt et al., 2001), pollutant removal (Ning et al., 2025; Osman et al., 2021), and other chemical systems (André et al., 2016). The solid materials in such systems, including various types of oxygen carriers (Cho et al., 2004) and adsorbents (Abanades et al., 2005; Lupiáñez et al., 2013), are commonly prepared as porous particles, which is especially applicable to fluidized bed reactors. Multiple physical and chemical processes simultaneously occur in the reactor, which leads to sophisticated impacts of design parameters on the reactor’s performance, requiring to be described by a whole process model.

A multiscale modeling framework has been proposed for catalytic reactions (Bruix et al., 2019), spanning from the microscopic atomic behaviors and surface microkinetics to the macroscopic flows; however, this framework does not cover non-catalytic reactions on the aspect of solid conversion, particularly for solids in internal lattices, compared with catalytic reactions during which the solid components remain unchanged. A complete framework for non-catalytic heterogeneous reactions can be constructed by integrating mesoscale mass transfer processes between microkinetics and fluidization, which consists of five scales as Figure 1, including:

1. The atom scale, where solid atoms and gas molecules interact, triggering an elementary reaction;

2. The surface scale, where reaction intermediates are adsorbed onto active sites distributed on the solid surface, according to Langmuir’s adsorption model (Langmuir, 1916);

3. The grain scale, where solid ions migrate between the surface and internal lattices of a grain, converting the solid reactant to the product;

4. The particle scale, where gas molecules diffuse through pores in a particle, before reaching the inner grain surfaces;

5. The reactor scale, where the particles are fluidized by the gas, forming a two-phase flow.

Figure 1
Diagram illustrating a scale from atom to reactor, showing an atom structure, molecular surface, grain composition, porous particle, and reactor filled with material. The scale indicates lengths and times from 10^-9 to 10^0 meters and seconds.

Figure 1. Multiscale modeling framework of non-catalytic heterogeneous reactions.

Various validated models for every independent scale are found in the literature. At the atom scale, the reaction paths and rate constants are directly calculated from the material properties with DFT (density functional theory), also known as the first-principles calculation (Bhandari et al., 2020; Fang et al., 2022; Li et al., 2023; Liu T. et al., 2021; Noor et al., 2020; Wang X. et al., 2024; Yan et al., 2020), which is a rigorous theoretical model requiring no empirical parameters. At the surface scale, the microkinetics, evolving the coverages of different species on the active sites, is modeled with mean-field approximation (Madon et al., 2011; Maestri et al., 2008; Thybaut et al., 2011) or kinetic Monte Carlo (Andersen et al., 2019; Reuter, 2016); elementary reactions may occur by single-site mechanisms including adsorption–desorption and Eley–Rideal, or multi-site mechanisms such as dissociation–association and Langmuir–Hinshelwood (Motagamwala and Dumesic, 2021); dual-site mechanisms are possible on a surface with different types of active sites (Van Belleghem et al., 2022; D’Ambrosio et al., 2024). At the grain scale, the solid conversion process can be described by homogeneous bulk diffusion (Arangio et al., 2015; Stearn and Eyring, 1940; Willis and Wilson, 2022), shrinking core model (Ishida and Wen, 1971; Szekely and Evans, 1970), or a more complex product island model (Fang et al., 2011; Li, 2020), based on the structure of the solid reactant. At the particle scale, intraparticle gas diffusion is controlled by the distribution of grains and pores, shaped by either grain models (Dam-Johansen et al., 1991; Gibson III and Harrison, 1980; Szekely and Evans, 1971) or pore models (Bhatia and Perlmutter, 1980; Bhatia and Perlmutter, 1981; He et al., 2013; Petersen, 1957; Sandmann Jr and Zygourakis, 1986). At the reactor scale, numerical simulation is applied to the two-phase flow to predict the fluidization behavior, the gas phase evolved with CFD (computational fluid dynamics), and the dense particle phase modeled with TFM (two-fluid model) (Ishii and Mishima, 1984), MP-PIC (multiphase particle-in-cell) (Andrews and O'Rourke, 1996) or DEM (discrete element method) (Cundall and Strack, 1979; Golshan et al., 2020); as the computational capacity develops, the most detailed and computationally expensive CFD–DEM model has been increasingly adopted, supporting up to 105–108 particles (Golshan et al., 2020).

Despite the variety of abovementioned models, the coupling between adjacent scales, which multiplies the computational costs of both models, has been an obstacle to completing the framework. Existing studies have proposed several coupling strategies for particular scales by simplifying the larger scale to identical sub-processes of the smaller scale. Microscopic coupling, suggested by studies on catalytic microkinetics (Alexopoulos et al., 2016; Chen and Wang, 2024; Jørgensen and Grönbeck, 2016; Rawal et al., 2021; Yin et al., 2020; Yu et al., 2023), assumes that active sites uniformly distribute on the surface, all sites sharing the same reaction paths derived from first-principles calculation; on this basis, the difference among sites is omitted in single-site reactions, while in multi-site reactions, the state of every pair of neighboring sites shall be considered (Razdan and Bhan, 2021). Macroscopic coupling appears in both CFD–DEM studies numerically solving intraparticle diffusion (Hadian et al., 2024), and particle-scale models analytically expressing the reaction rate with the surface gas concentration (Sedghkerdar and Mahinpey, 2015; Wang et al., 2017; Wang et al., 2021; Yang et al., 2016). The above micro- and macroscopic coupling strategies have been generally agreed in the literature; in contrast, mesoscopic coupling, focusing on the relation between microkinetics and solid conversion, has not undergone sufficient research. Some studies have introduced a bulk-phase ion diffusion process, associating the surface coverage of active ions with its bulk concentration (Li, 2022; Li L. et al., 2021; Liu et al., 2022; Wang et al., 2025; Wang Y. et al., 2024); others have described a phase transformation process based on the shrinking core model, which is integrated into the microkinetics (Cai and Li, 2024a; Cai and Li, 2024b; Cai and Li, 2024c).

However, none of above coupling models have accomplished a complete multiscale framework as Figure 1. Studies on microkinetics have not considered the solid conversion processes of a real porous particle; CFD–DEM studies, mainly applied to the combustion of organic solid fuels, whose molecular structures are unclear, have not developed verifiable microkinetics. Consequently, parameters from the absent scales rely on experimental fitting, instead of theoretical calculation or direct measurement. Models with mesoscopic coupling, excluding the reactor-scale, are unable to provide adequate information of the diverse particles; furthermore, only surfaces with a single type of active sites are involved in these studies (Cai and Li, 2024a; Cai and Li, 2024b; Cai and Li, 2024c; Li, 2022; Li T. et al., 2021; Liu et al., 2022; Wang et al., 2025; Wang Y. et al., 2024), while dual-site reaction mechanisms have not been discussed, resulting in a limited scope of application.

The aim of this study is to develop a multiscale model for non-catalytic heterogeneous reactions under the framework of Figure 1, coupling strategies introduced between every two adjacent scales. The atom-scale reaction path is derived from first-principles calculation; the surface-scale mean-field microkinetics describes the coverages of two types of active sites; the grain-scale ion diffusion is coupled with the microkinetics by the lumped parameter method; the particle-scale gas diffusion is treated under a uniform grain model; the reactor-scale fluidization is simulated with CFD–DEM. Experimental validation is conducted on a micro-fluidized-bed thermogravimetric analyzer (MFB–TGA) (Li et al., 2019). Without loss of generality, the dual-site reaction of a perovskite oxygen carrier (CaMn0.375Ti0.5Fe0.125O3−δ, CMTF8341) reduced by H2 is considered. Kinetics of this reaction has been experimentally studied in a prior work (Liu L. et al., 2021), where the material is prepared by spray drying and calcination, sieved between 180 and 250 μm, and tested on MFB–TGA; fast reaction kinetics and good stability have been observed under H2/O2 redox cycles. The gas reactant, H2, offers a special dual-site mechanism as needed, while fewer elementary reactions are involved compared with other reducing agents. Moreover, the perovskite structure allows oxygen anions to undergo homogeneous bulk diffusion, rather than changing the surface structure into another phase. Thus, this reaction is expected to concisely and clearly describe the full modeling framework.

2 Model

The overall reaction of H2 reducing CMTF8341 as Equation 1.

H2+1δCaMn0.375Ti0.5Fe0.125O3H2O+1δCaMn0.375Ti0.5Fe0.125O3δ(1)

Complete reduction is achieved when δ=0.5, corresponding to a solid mass loss of 5.73%; however, the actual mass loss, as measured in Section 3.1, is less than the theoretical value due to the impurity in CMTF8341. Given that the solid mass decreases from mox to a limit of mre, when the reactive component is completely reduced, the mass capacity of the oxygen carrier material is defined as Equation 2.

ROC=moxmremox(2)

The solid conversion is defined by the real-time solid mass, mt, as Equation 3.

X=moxmtmoxROC(3)

which is applicable to a grain, a particle, or all CMTF8341 particles in a reactor. The conversion increases from 0 to 1 during the whole reaction process. The final output of the model is the change of conversion against time in a reactor, which is also measured by the experiment as Section 3.

2.1 Atom scale

DFT calculation is conducted for H2 reducing fully oxidized CMTF8341. The reaction path is illustrated in Figure 2, consisting of (three elementary steps) as Equations 46. First, one H2 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

H2g+O*+#OH*+H#(4)

Figure 2
Atomic diagrams depict the reduction process on a CMTF8341 lattice surface, consisting of Ca, Mn, Ti, Fe, O and H atoms. Active oxygen is highlighted as green. Reaction occurs in the sequence of: (1) H₂+O*+#, (2) OH*+H#, (3) H₂O*+#, and (4) H₂O+*+#.

Figure 2. Reaction path of CMTF8341 reduced by H2 from DFT calculation. The H2 molecule undergoes asymmetrical dissociation, re-association, and desorption to become H2O.

Subsequently, the single H atom is associated with the OH radical, forming one H2O molecule on the Mn site.

OH*+H#H2O*+#(5)

The last step is the H2O molecule desorbing into the gas phase, leaving a vacant Mn site on the surface.

H2O*H2Og+*(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.

Figure 3
Reaction coordinate diagram showing the relative energy changes during a chemical process. Initial state H₂+O*+# (0 eV). Transition state TS1 (1.01 eV). Intermediate state OH*+H# (-2.82 eV). Transition state TS2 (-2.79 eV). Intermediate state H₂O*+# (-3.32 eV). Final state H₂O+*+# (-2.52 eV).

Figure 3. Energy diagram of CMTF8341 reduced by H2 from DFT calculation.

Table 1
www.frontiersin.org

Table 1. Properties of species and transition states involved in the reaction path.

2.2 Surface scale

2.2.1 Mean-field assumption

The surface of the CMTF8341 crystal is highly periodical, where both Mn sites and Ca sites are evenly distributed. A mean-field assumption is employed in treating the elementary reactions on different sites, which neglects the change of the surface force field caused by different adsorbates, so that all sites are in an identical environment. Consequently, elementary reactions on every site occur along the same paths as Figure 2, while every site shares an equal reacting probability.

Previous studies on single-site reactions, based on the mean-field assumption, have treated the surface species as non-localized independent particle systems (Cai and Li, 2024a; Cai and Li, 2024b; Cai and Li, 2024c; Liu et al., 2022; Wang et al., 2025). A dual-site reaction, however, can only occur on an adjacent pair of sites, rather than every pair; thus, the surface species shall not be treated as completely independent for dual-site reactions. Alternatively, a pair of sites can be regarded independent of other pairs, given that one elementary reaction does not interfere with reactions on other sites.

Consider a CMTF8341 surface, where the numbers of Mn sites and Ca sites are N*,tot and N#,tot, respectively. Every Mn site has an equal number of neighboring Ca sites, the number denoted as d*, while every Ca site has d# neighboring Mn sites. The total number of Mn–Ca site pairs on the surface is thus as Equation 7.

N*#,tot=N*,totd*=N#,totd#(7)

Assuming that a Mn site is occupied by a surface species X*, and its neighboring Ca site by Y#, the state of this site pair can be expressed as X*Y#. As X* and Y# randomly distribute on their corresponding sites, the number of X*Y# pairs has a mathematical expectation of Equation 8

NX*Y#=N*#,totθX*ϑY#(8)

where θX* is the coverage of X* on Mn sites, and ϑY# that of Y# on Ca sites.

2.2.2 Rate equation

The rate of an elementary reaction is derived from the transition state theory. The transition state of the single-site reaction (Equation 6) is formed on a Mn site, so the rate equation is based on the Mn-site coverages, as Equation 9.

θ˙ER3=N˙ER3N*,tot=k+3θH2O*k3pH2Oθ*k+3=kBThqTS30,qH2O*0expΔε0,+3kBTk3=1hqTS30,qH2O0Vq*0expΔε0,3kBT(9)

where every q0 represents the partition function of a species (or transition state) derived from statistical mechanics; the energy barriers, Δε0, are given by Figure 3.

As for a dual-site reaction, taking Equation 4 for example, the transition state is formed on a site pair, so the rate is expressed on a site-pair basis, as Equation 10

N˙ER1N*#,tot=k+1pH2θO*ϑ#k1θOH*ϑH#(10)

where N˙ER1 is the net turnover frequency of Equation 4 on the entire surface. Whereas the solid reactant and product species (O* and *) are located on Mn sites, the rate equation is further expressed with the Mn-site coverage. Given Equation 11

θ˙ER1=N˙ER1N*,tot(11)

the rate equation is transformed into Equation 12.

θ˙ER1=d*k+1pH2θO*ϑ#d*k1θOH*ϑH#k+1=1hqTS10,qH20VqO*0q#0expΔε0,+1kBTk1=kBThqTS10,qOH*0qH#0expΔε0,1kBT(12)

The above derivation is the same for Equation 5, resulting in Equation 13.

θ˙ER2=d*k+2θOH*ϑH#d*k2θH2O*ϑ#k+2=kBThqTS20,qOH*0qH#0expΔε0,+2kBTk2=kBThqTS20,qH2O*0q#0expΔε0,2kBT(13)

In brief, the rate constants of dual-site reactions shall be multiplied by the number of neighbors (d*) when considering one type of sites. This result corresponds to the theory by Nørskov et al. (2014), which suggests that every neighbor provides an equivalent reaction path, thus increasing the reacting probability by a factor of d*.

2.3 Grain scale

2.3.1 Bulk ion diffusion

The reduction of CMTF8341 is achieved by the removal of oxygen anions (O2) in the lattices, while only surface oxygen (O*) is consumed through Equation 4. Lattice oxygen does not directly participate in the surface reactions; instead, it diffuses outward from the internal lattices to the grain surface, driven by its concentration gradient. The structure of CMTF8341 remains stable during the process, forming a solid solution where lattice oxygen (OO) and vacancies (VO) act as crystal defects, rather than mixed phases of crystals.

The ion diffusion process is fast enough in perovskites (Li, 2022), so that the concentration gradient is negligible. A lumped parameter method is employed by treating OO and VO as two non-localized independent particle systems. The surface–lattice oxygen transformation can be expressed as

OO+*VO+O*(14)

The equilibrium of Equation 14 is reached by equalizing the surface and lattice oxygen concentrations, given by Equation 15

θO*θ*=COOCVO(15)

where C is the concentration of OO or VO. The diffusion rate is expressed as Equation 16

θ˙ER4=JΛ(16)

where J is the surface diffusion flux, and Λ the areal Mn-site density.

The overall reaction is the sum of elementary reactions Equations 4, 5, 6, 14, as Equation 17

H2g+OOH2Og+VO(17)

which is equivalent to Equation 1.

2.3.2 Non-catalytic microkinetics

Given the rates of surface reactions Equations 46 and ion diffusion Equation 14, the complete microkinetics on Mn sites is expressed as Equation 18

dθO*dt=θ˙ER1+θ˙ER4dθOH*dt=θ˙ER1θ˙ER2dθH2O*dt=θ˙ER2θ˙ER3dθ*dt=θ˙ER3θ˙ER4(18)

while on Ca sites, there is Equation 19.

dϑH#dt=N*,totN#,totθ˙ER1θ˙ER2dϑ#dt=N*,totN#,totθ˙ER1+θ˙ER2(19)

For lattice oxygen or vacancies, there is Equation 20.

43πrg3dCOOdt=4πrg2Λθ˙ER443πrg3dCVOdt=4πrg2Λθ˙ER4(20)

Note that the ion diffusion rate, θ˙ER4, appears only in Equations 18, 20, which indicates that bulk ion diffusion occurs only on Mn sites, rather than Ca sites. Therefore, Mn is the primary site in this non-catalytic microkinetics, directly affecting the conversion of a grain; Ca is the secondary site accommodating intermediates which assist the reaction process.

The above equations result in a total of eight species evolved by four elementary reactions. The number of independent variables among the species shall be no more than the number of reactions. One restriction is that the sum of all quantities on one type of sites (or the lattices) shall be constant, as Equation 21

θO*+θOH*+θH2O*+θ*=1ϑH#+ϑ#=1COO+CVO=Cmax=ρtROCMO(21)

where ρt is the solid true density, and MO is the molar mass of O2. Meanwhile, OH* and H# are always formed or consumed in pairs, giving Equation 22.

ϑH#=N*,totN#,totθOH*(22)

Consequently, the state of a grain can be determined by four independent variables (for example, θO*, θOH*, θH2O*, and COO), which equals the number of elementary reactions.

2.3.3 Quasi-steady approximation

Further simplification is obtained from the partial equilibrium assumption, suggesting that the overall reaction rate is controlled by the slowest step, which is Equation 4 in this case; other reversible elementary reactions, Equations 5, 6, 14, are near equilibrium, as concluded by the results in Section 4.2. However, no analytic solution is obtained from the above assumption. Given that the intermediates, OH* and H2O*, undergo rapid generation and consumption, resulting in Equation 23

θOH*1θH2O*1(23)

an approximate solution is derived as Equation 24

θ*1COOCmaxθO*COOCmaxθH2O*1K3pH2O1COOCmaxθOH*1K2K3pH2O1COOCmax(24)

where Ki=k+i/ki is the equilibrium constant.

As only one independent variable (COO) is left among all species, the conversion of a grain is expressed as Equation 25.

Xg=1COOCmax(25)

Substitute Xg into Equations 18, 20, eliminate the rates of Equations 5, 6, 14, and approximate the rate of Equation 4 with its forward rate, giving

dXgdt11+ρtROCrg3MOΛ·d*k+11XgpH2(26)

The preceding factor of Equation 26 indicates the effect of ion diffusion compared with surface reaction, hereinafter expressed with the symbol Li proposed by Li Z. et al. (2021) and Li (2022), as Equation 27.

Li=11+ρtROCrg3MOΛ(27)

2.4 Particle scale

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

pt=1r2rr2Deffprωpprr=0=0pr=rp=ps(28)

where rp is the particle radius, and ps is the surface gas partial pressure. The reaction term is derived from the gas–solid stoichiometry, as Equation 29.

ω=1ααρtROCRTMO·Li·k+11Xp(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.

dXpdt=0rpLi·k+11Xgp·4πr2dr/43πrp3(30)

The governing equation, Equation 28, requires numerical solution due to both p and Xg 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

dXpdt=η·Li·k+11Xpps(31)

where η is the effectiveness factor of intraparticle gas diffusion, given by Equation 32.

η=3ϕ2ϕcothϕ1ϕ=rpωDeff(32)

the Thiele modulus ϕ representing the influence of reaction versus diffusion.

2.5 Reactor scale

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.

2.5.1 Particle phase

The particle mass is related to its conversion by Equation 3; consequently, the mass equation of a particle is Equation 33

dmpdt=mp,oxROCdXpdt(33)

where mp,ox is the mass of a fully oxidized particle, and the conversion rate is given by Equation 31.

The translational motion of a particle is evolved by Newton’s Second Law as Equations 34, 35

dxpdt=vp(34)
mpdvpdt=fdrag+fcoll+mpg(35)

and rotational motion as Equation 36

Ipdwpdt=Tcoll(36)

where xp is the position, vp the velocity, and wp the angular velocity. The drag force, fdrag, is usually expressed in the form of Equation 37

fdrag=π6dp3βpufvp(37)

where ufvp is the fluid velocity relative to the particle, and dp is the particle diameter; the coefficient βp is given by the Koch–Hill model, applicable to both dilute and dense phases (Koch and Hill, 2001). The collision force, fcoll, is decomposed into a normal and a tangential component by the spring-slider-dashpot model (Cundall and Strack, 1979), as Equation 38

fN=κNΔxNγNΔvNfT=minκTΔxTγTΔvT,μfN(38)

where the subscripts “N” and “T” represent the normal and tangential components, respectively. Both components consist of an elastic force proportional to the overlap Δx, and a damping force depending on the relative velocity Δv, while fT shall not exceed the sliding friction μfN. The elastic and damping coefficients, κ and γ, are derived from the restitution coefficient, Young’s modulus, and Poisson’s ratio of the material. Detailed expressions can be found in the computational software (Kloss et al., 2012).

2.5.2 Fluid phase

The gas flow in the reactor is regarded incompressible, described by the governing equations of fluid dynamics. As part of the volume is occupied by the dense particle flow, the governing equations shall involve the fluid volume fraction, denoted as αf. Consequently, the equations of continuity and species transfer are given by Equation 39

αft+·αfuf=1ρfVcelldmpdtαfρjt+·αfufρj·αfDjρj=1Vcellm˙p,j(39)

where ρj is the partial density of gas species j. Interphase mass transfer is introduced by the extra source terms, containing the particle reaction rates as Equation 33.

The equation of momentum is expressed as Equation 40

αfuft+·αfufuf·τ=pρf+αfg1ρfVcellfdrag+dmpdtufτ=νeffuf+ufT23·ufI(40)

where the interphase drag force given by Equation 37 is considered, along with the momentum carried by the transferred mass. The viscosity νeff contains both laminar and turbulent components, with turbulence solved by the kϵ model.

3 Experimental and computational setup

3.1 Experiment

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 (N2), oxidizing (O2) and reducing (H2) components, under different H2 concentrations (5 vol%, 10 vol% and 20 vol%) and temperatures (750 °C, 800 °C, 850 °C and 900 °C).

Table 2
www.frontiersin.org

Table 2. Fluidization properties of particles.

3.2 Computation

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.

Table 3
www.frontiersin.org

Table 3. Parameters used in CFD–DEM simulation.

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.

Figure 4
Class diagram illustrating relationships and attributes for chemical reaction modeling. It includes classes like ElementaryReaction, TransitionState, StableSpecies, and FirstPrinciplesRateEquation, with attributes and methods related to energy, state, and calculation of reaction rates and heats. Inheritance and associations are depicted, indicating hierarchical structure and shared attributes across entities such as Species, State, Translation, Vibration, and Rotation classes.

Figure 4. Architecture of the single-particle reaction kinetics program. The abstract class Species is used to construct different species and transition states in the reaction path; statistical-mechanical calculation methods are implemented; elementary and overall reactions are constructed from given species.

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.

4 Results and discussion

4.1 Solid conversion rate

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.

Figure 5
Two graphs labeled (A) and (B) show the relationship between conversion (X) and time (t/s). Graph (A) is for 900°C, depicting experimental and simulated data for 5, 10, and 20 volume percent. Graph (B) is for 10 volume percent, comparing experimental and simulated data at 750°C, 800°C, 850°C, and 900°C. Both graphs indicate various curves and data points to highlight experimental and simulation results.

Figure 5. Solid conversion results from computation and experiment at (A) different H2 concentrations, (B) different temperatures.

The conversion curves under different inlet H2 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% H2, while the time is prolonged to ∼15 s for 10 vol% H2, and ∼30 s for 5 vol% H2. 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 H2 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 H2 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.

Figure 6
Logarithmic graph depicting the relationship of rate constant, exponential factor, and partition function quotient versus inverse temperature. The black line represents the rate constant, the red line represents the exponential factor, and the blue line represents the partition function quotient. Temperature is marked in Kelvin along the X-axis from zero to five times ten to the power of negative three inverse Kelvin. The Y-axis scale ranges from ten to the power of negative forty to ten to the power of ten in scientific notation. The graph includes temperature markers at one thousand, five hundred, and two hundred fifty Kelvin.

Figure 6. Change of the rate constant and its factors (in logarithm) against inversed 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 H2 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).

Figure 7
Scatter plot comparing experimental (X(Exp)) and simulated (X(Sim)) values, ranging from zero to one. Symbols represent different conditions like volume percentage and temperature, with a diagonal reference line and two dashed lines. Legend identifies specific conditions using distinct symbols.

Figure 7. Error margin of computed solid conversion compared with experimental results.

4.2 Microkinetics with bulk diffusion

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−8 mol/m2) 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 1820 is originally solved as an ordinary differential equation set; an example condition of 20 vol% H2, 1 vol% H2O 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 H2O*) 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*, H2O*, and * in order. An equilibrium is established for reaction Equation 6 immediately after the reaction begins, owing to a high turnover frequency of the H2O 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.

Figure 8
Three graphs are displayed: (A) shows the concentration change of O* (red) decreasing and * (black) increasing over 10 seconds. (B) illustrates the increasing trend of OH* (red) and H2O* (blue) concentrations over 10 seconds. (C) depicts rate changes for three reactions, ER1, ER2, and ER3, with forward and reversed rates, showing different levels of stabilization over 10 seconds.

Figure 8. Numerical solution to the dual-site microkinetics with bulk diffusion. (A) Coverages of reactant and product species; (B) Coverages of intermediates; (C) Forward and reversed elementary reaction rates.

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.

pH2OCVOpH2COO=qH2O0Vq*0qH20VqO*0expΔε0,overallkBT(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 pH2 and pH2O have an order of kPa–MPa. Only under extreme conditions when pH2/pH2O<1013, typically when the gas reactant is depleted due to insufficiency, can the solid achieve incomplete conversion.

Figure 9
Graph showing the equilibrium fraction \(X_{g,eq}\) versus the ratio \(p_{H2}/p_{H2O}\) for various temperatures: 750°C, 800°C, 850°C, and 900°C. Each temperature is represented by a different colored curve with a steep increase around \(10^{-15}\) on the x-axis.

Figure 9. Change of equilibrium solid conversion against gas concentrations.

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% H2 and 900 °C, as d* increases from 1 to 4. Rates proportional to d* are observed in the early stage when Xg<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.

Figure 10
Graph showing four curves representing \(X_g\) over time \(t/s\) for different values of \(d_*\). Blue curve for \(d_* = 1\), red for \(d_* = 2\), yellow for \(d_* = 3\), and purple for \(d_* = 4\). All curves increase, with the purple curve reaching the highest value.

Figure 10. Effect of neighboring-site number on the grain’s conversion rate (20 vol% H2, 900 °C).

4.3 Particle diversity

The evolution of particle distribution at 900 °C, 20 vol% H2, 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.

Figure 11
Seven sequential frames depict granular mixing in a cylindrical container from 0 to 12 seconds. Red and blue particles mix over time, starting densely packed and gradually spreading with increasing uniformity. The color gradient indicates particle concentration.

Figure 11. Particle location and conversion distributions at different moments.

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% H2, 900 °C. Particle 3 undergoes an average-rate reaction, while Particle 1 is faster and Particle 2 slower. These particles are exposed to particular H2 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 H2 than Particles 3 and 2.

Figure 12
Panel (A) shows a line graph with conversion versus time in seconds. It includes curves for average, maximum, and minimum conversion, along with data points for three particles. Panel (B) is a three-dimensional scatter plot in a cylinder, displaying the positions of three different particles, indicated by blue, red, and green dots. Panel (C) repeats the same plot with emphasis on cross-sections for each particle group.

Figure 12. Conversion of three typical particles undergoing fast (particle 1), slow (particle 2) and average (particle 3) reactions. (A) Conversion curves; (B) Positions at different moments; (C) Average positions during the initial 5 s.

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 H2 is consumed by the few particles near the inlet, while other particles only obtain a low H2 concentration.

Figure 13
(A) Six histograms showing the distribution of a variable \( X \) over time, from \( t = 1.0 \) to \( 9.0 \) seconds, with increasing average values. (B) A line graph depicting standard deviation versus average conversion for different conditions, with multiple lines representing different volume percentages and temperatures. (C) A scatter plot showing maximum standard deviation against theoretical rate, with a linear trend line.

Figure 13. Standard deviation of particle conversion at different moments. (A) Histograms; (B) Standard deviation against average conversion; (C) Maximum standard deviation against theoretical reaction rate.

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 H2 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.

4.4 Intraparticle gas diffusion

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 H2 has a relatively low rate compared with gas diffusion, so intraparticle gas diffusion has a minor impact among all processes in the multiscale framework.

Figure 14
Line graph showing efficiency (η) versus conversion (Xₚ) for temperatures 750 °C, 800 °C, 850 °C, and 900 °C. Efficiency increases with conversion for all temperatures. Higher temperatures show slightly lower efficiency increases, with 750 °C being the highest and 900 °C the lowest.

Figure 14. Effectiveness factor of intraparticle gas diffusion changing with particle conversion.

4.5 Overall reaction order

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 1Xg. 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 H2 reduction kinetics (Li and Li, 2024; Li, 2022; Li Z. et al., 2021) have generally adopted the symmetrical dissociation mechanism as Equations 4244

H2g+2O*2OH*(42)
2OH*H2O*+O*(43)
H2O*H2Og+*(44)

where the H2 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

dXgdt=Li·d*k+11Xg2pH2(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.

4.6 Computational cost for scale-up systems

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.

Figure 15
Line graph showing the relationship between particles per CPU and real-time per simulation. The x-axis represents particles per CPU on a logarithmic scale, while the y-axis represents real-time per second of simulation. Data points are labeled with numbers indicating CPU counts: 48, 32, 24, 12, 8, 4, and 2. A minimum real-time appears at 32 CPU cores.

Figure 15. Computational time varying against the number of particles per processor. The processor number is tagged beside each data point.

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.

5 Conclusion

A multiscale model is developed regarding the whole process of non-catalytic heterogeneous reactions in a fluidized bed. Physical and chemical processes spanning across various scales are involved, including the atom-scale elementary reaction path, the surface-scale dual-site microkinetics, the grain-scale asymmetrical bulk diffusion, the particle-scale gas diffusion, and the reactor-scale fluidization behavior. The model is applied to the reduction of CMTF8341 by H2 and validated by the MFB–TGA experiment, revealing the effects of inlet gas concentration and temperature on the macroscopic reaction kinetics.

The prediction of overall reaction rate strongly relies on the reaction path from DFT calculation, whose results suggest a dual-site dissociation for H2 molecules, resulting in a first-order overall rate equation compared with the second-order symmetrical dissociation mechanism. A mean-field microkinetics regarding active site pairs is established, every site connected to an equal number of neighbors on the periodical surface, forming the same number of equivalent reaction paths. Furthermore, the dual-site microkinetics is coupled with bulk diffusion in an asymmetrical approach, where the primary site exchanges ions with inner lattices, while the secondary site only accommodates intermediates to assist surface reactions. The neighboring-site number is currently a hypothetic parameter; scanning transmission electron microscopic (STEM) and scanning tunneling microscopic (STM) observation are potential ways to justify corresponding assumptions.

The rigorous theoretical model is based on physical parameters instead of experimental fitting, thus reducing the experimental costs when applied to other reaction systems. An extensible software package is developed for mesoscopic models and integrated into open source CFD–DEM software; generalization to other heterogeneous reactions can be accomplished by inputting DFT results along with elementary reaction equations and instrumentally measured properties. This work is conducted on a lab-scale reactor, while the computational cost is an obstacle to scaling the model to industrial-scale systems; potential methods including coarse-graining models and GPU acceleration may be involved in relevant future studies.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

RW: Data curation, Methodology, Formal Analysis, Investigation, Software, Conceptualization, Writing – original draft. ZL: Writing – review and editing, Funding acquisition, Conceptualization, Supervision. LL: Validation, Methodology, Writing – review and editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported by the National Key Research and Development Program of China (No. 2024YFB4105904).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2025.1656180/full#supplementary-material

References

Abanades, J. C., Anthony, E. J., Wang, J., and Oakey, J. E. (2005). Fluidized bed combustion systems integrating CO2 capture with CaO. Environ. Sci. Technol. 39 (8), 2861–2866. doi:10.1021/es0496221

PubMed Abstract | CrossRef Full Text | Google Scholar

Alexopoulos, K., John, M., Van der Borght, K., Galvita, V., Reyniers, M. F., and Marin, G. B. (2016). DFT-based microkinetic modeling of ethanol dehydration in H-ZSM-5. J. Catal. 339, 173–185. doi:10.1016/j.jcat.2016.04.020

CrossRef Full Text | Google Scholar

Andersen, M., Panosetti, C., and Reuter, K. (2019). A practical guide to surface kinetic Monte Carlo simulations. Front. Chem. 7, 202. doi:10.3389/fchem.2019.00202

PubMed Abstract | CrossRef Full Text | Google Scholar

André, L., Abanades, S., and Flamant, G. (2016). Screening of thermochemical systems based on solid-gas reversible reactions for high temperature solar thermal energy storage. Renew. Sust. Energy Rev. 64, 703–715. doi:10.1016/j.rser.2016.06.043

CrossRef Full Text | Google Scholar

Andrews, M. J., and O'Rourke, P. J. (1996). The multiphase particle-in-cell (MP-PIC) method for dense particulate flows. Int. J. Multiph. Flow. 22 (2), 379–402. doi:10.1016/0301-9322(95)00072-0

CrossRef Full Text | Google Scholar

Arangio, A. M., Slade, J. H., Berkemeier, T., Pöschl, U., Knopf, D. A., and Shiraiwa, M. (2015). Multiphase chemical kinetics of OH radical uptake by molecular organic markers of biomass burning aerosols: humidity and temperature dependence, surface reaction, and bulk diffusion. J. Phys. Chem. A 119 (19), 4533–4544. doi:10.1021/jp510489z

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhandari, S., Rangarajan, S., Maravelias, C. T., Dumesic, J. A., and Mavrikakis, M. (2020). Reaction mechanism of vapor-phase formic acid decomposition over platinum catalysts: DFT, reaction kinetics experiments, and microkinetic modeling. ACS Catal. 10 (7), 4112–4126. doi:10.1021/acscatal.9b05424

CrossRef Full Text | Google Scholar

Bhatia, S. K., and Perlmutter, D. D. (1980). A random pore model for fluid-solid reactions: I. isothermal, kinetic control. AIChE J. 26, 379–386. doi:10.1002/aic.690260308

CrossRef Full Text | Google Scholar

Bhatia, S. K., and Perlmutter, D. D. (1981). A random pore model for fluid-solid reactions: II. diffusion and transport effects. AIChE J. 27, 247–254. doi:10.1002/aic.690270211

CrossRef Full Text | Google Scholar

Bruix, A., Margraf, J. T., Andersen, M., and Reuter, K. (2019). First-principles-based multiscale modelling of heterogeneous catalysis. Nat. Catal. 2 (8), 659–670. doi:10.1038/s41929-019-0298-3

CrossRef Full Text | Google Scholar

Cai, J., and Li, Z. (2024a). First principle-based rate equation theory for the carbonation kinetics of CaO with CO2 in calcium looping. Chem. Eng. J. 479, 147484. doi:10.1016/j.cej.2023.147484

CrossRef Full Text | Google Scholar

Cai, J., and Li, Z. (2024b). First principle-based rate equation theory for the calcination kinetics of CaCO3 in calcium looping. Energy fuels. 38 (9), 8145–8156. doi:10.1021/acs.energyfuels.4c00877

CrossRef Full Text | Google Scholar

Cai, J., and Li, Z. (2024c). First principles-based kinetic analysis of Ca(OH)2 dehydration in thermochemical energy storage system. J. Energy Storage 97, 112759. doi:10.1016/j.est.2024.112759

CrossRef Full Text | Google Scholar

Chen, Y. C., and Wang, G. C. (2024). Structure sensitivity of the water–gas shift reaction over Cu/CeO2: a combination of DFT, MF-MKM, and kMC. J. Phys. Chem. C 128 (24), 9926–9939. doi:10.1021/acs.jpcc.4c01539

CrossRef Full Text | Google Scholar

Cho, P., Mattisson, T., and Lyngfelt, A. (2004). Comparison of iron-nickel-copper-and manganese-based oxygen carriers for chemical-looping combustion. Fuel 83 (9), 1215–1225. doi:10.1016/j.fuel.2003.11.013

CrossRef Full Text | Google Scholar

Christensen, R., Hansen, H. A., and Vegge, T. (2015). Identifying systematic DFT errors in catalytic reactions. Catal. Sci. Technol. 5, 4946–4949. doi:10.1039/C5CY01332A

CrossRef Full Text | Google Scholar

Cundall, P. A., and Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Géotechnique 29 (1), 47–65. doi:10.1680/geot.1979.29.1.47

CrossRef Full Text | Google Scholar

Dam-Johansen, K., Hansen, P. F. B., and Østergaard, K. (1991). High-temperature reaction between sulphur dioxide and limestone—III. A grain-micrograin model and its verification. Chem. Eng. Sci. 46 (3), 847–853. doi:10.1016/0009-2509(91)80191-Z

CrossRef Full Text | Google Scholar

Derksen, J. J. (2014). Simulations of solid–liquid mass transfer in fixed and fluidized beds. Chem. Eng. J. 255, 233–244. doi:10.1016/j.cej.2014.06.067

CrossRef Full Text | Google Scholar

D’Ambrosio, A., Bertino, A., Todaro, S., Santoro, M., Cannilla, C., Frusteri, F., et al. (2024). Kinetic modeling of the direct dimethyl ether (DME) synthesis over hybrid multi-site catalysts. Catalysts 14 (1), 61. doi:10.3390/catal14010061

CrossRef Full Text | Google Scholar

Fang, F., Li, Z. S., Cai, N. S., Tang, X. Y., and Yang, H. T. (2011). AFM investigation of solid product layers of MgSO4 generated on MgO surfaces for the reaction of MgO with SO2 and O2. Chem. Eng. Sci. 66 (6), 1142–1149. doi:10.1016/j.ces.2010.12.014

CrossRef Full Text | Google Scholar

Fang, D., Hou, S., Ye, Y., Jin, Q., He, F., and Xie, J. (2022). Insight into highly efficient FeOx catalysts for the selective catalytic reduction of NOx by NH3: experimental and DFT study. Appl. Surf. Sci. 599, 153998. doi:10.1016/j.apsusc.2022.153998

CrossRef Full Text | Google Scholar

García-Labiano, F., de Diego, L. F., Adánez, J., Abad, A., and Gayán, P. (2005). Temperature variations in the oxygen carrier particles during their reduction and oxidation in a chemical-looping combustion system. Chem. Eng. Sci. 60, 851–862. doi:10.1016/j.ces.2004.09.049

CrossRef Full Text | Google Scholar

Gibson III, J. B., and Harrison, D. P. (1980). The reaction between hydrogen sulfide and spherical pellets of zinc oxide. Ind. Eng. Chem. Process Des. Dev. 19 (2), 231–237. doi:10.1021/i260074a005

CrossRef Full Text | Google Scholar

Golshan, S., Sotudeh-Gharebagh, R., Zarghami, R., Mostoufi, N., Blais, B., and Kuipers, J. A. M. (2020). Review and implementation of CFD-DEM applied to chemical process systems. Chem. Eng. Sci. 221, 115646. doi:10.1016/j.ces.2020.115646

CrossRef Full Text | Google Scholar

Hadian, M., de Munck, M. J. A., Buist, K. A., Bos, A. N. R., and Kuipers, J. A. M. (2024). Modeling of a catalytic fluidized bed reactor via coupled CFD-DEM with MGM: from intra-particle scale to reactor scale. Chem. Eng. Sci. 284, 119473. doi:10.1016/j.ces.2023.119473

CrossRef Full Text | Google Scholar

Hallberg, P., Leion, H., and Lyngfelt, A. (2011). A method for determination of reaction enthalpy of oxygen carriers for chemical looping combustion – application to ilmenite. Thermochim. Acta 524, 62–67. doi:10.1016/j.tca.2011.06.015

CrossRef Full Text | Google Scholar

He, W., Liu, Y., He, R., Ito, T., Suda, T., Fujimori, T., et al. (2013). Combustion rate for char with fractal pore characteristics. Combust. Sci. Technol. 185 (11), 1624–1643. doi:10.1080/00102202.2013.822370

CrossRef Full Text | Google Scholar

Ishida, M., and Jin, H. (1994). A new advanced power-generation system using chemical-looping combustion. Energy 19 (4), 415–422. doi:10.1016/0360-5442(94)90120-1

CrossRef Full Text | Google Scholar

Ishida, M., and Wen, C. Y. (1971). Comparison of zone-reaction model and unreacted-core shrinking model in solid–gas reactions—I Isothermal analysis. Chem. Eng. Sci. 26 (7), 1031–1041. doi:10.1016/0009-2509(71)80017-9

CrossRef Full Text | Google Scholar

Ishii, M., and Mishima, K. (1984). Two-fluid model and hydrodynamic constitutive relations. Nucl. Eng. Des. 82 (2–3), 107–126. doi:10.1016/0029-5493(84)90207-3

CrossRef Full Text | Google Scholar

Jørgensen, M., and Grönbeck, H. (2016). First-principles microkinetic modeling of methane oxidation over Pd(100) and Pd(111). ACS Catal. 6 (10), 6730–6738. doi:10.1021/acscatal.6b01752

CrossRef Full Text | Google Scholar

Kim, M. C., Sim, E., and Burke, K. (2013). Understanding and reducing errors in density functional calculations. Phys. Rev. Lett. 111, 073003. doi:10.1103/PhysRevLett.111.073003

PubMed Abstract | CrossRef Full Text | Google Scholar

Kloss, C., Goniva, C., Hager, A., Amberger, S., and Pirker, S. (2012). Models, algorithms and validation for opensource DEM and CFD-DEM. Prog. Comput. Fluid Dy. 12 (2–3), 140–152. doi:10.1504/pcfd.2012.047457

CrossRef Full Text | Google Scholar

Koch, D. L., and Hill, R. J. (2001). Inertial effects in suspension and porous-media flows. Annu. Rev. Fluid Mech. 33, 619–647. doi:10.1146/annurev.fluid.33.1.619

CrossRef Full Text | Google Scholar

Kresse, G., and Furthmüller, J. (1996). Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 54, 11169–11186. doi:10.1103/PhysRevB.54.11169

PubMed Abstract | CrossRef Full Text | Google Scholar

Langmuir, I. (1916). The constitution and fundamental properties of solids and liquids. Part I. Solids. J. Am. Chem. Soc. 38 (11), 2221–2295. doi:10.1021/ja02268a002

CrossRef Full Text | Google Scholar

Li, Z. (2020). General rate equation theory for gas–solid reaction kinetics and its application to CaO carbonation. Chem. Eng. Sci. 227, 115902. doi:10.1016/j.ces.2020.115902

CrossRef Full Text | Google Scholar

Li, Z. (2022). First-principles-based microkinetic rate equation theory for oxygen carrier reduction in chemical looping. Chem. Eng. Sci. 247, 117042. doi:10.1016/j.ces.2021.117042

CrossRef Full Text | Google Scholar

Li, J., and Li, Z. (2024). First principle based rate equation (1pRE) for reduction kinetics of Fe2O3 with syngas in chemical looping. Proc. Combust. Inst. 40 (1–4), 105363. doi:10.1016/j.proci.2024.105363

CrossRef Full Text | Google Scholar

Li, Y., Wang, H., Li, W., Li, Z., and Cai, N. (2019). CO2 gasification of a lignite char in microfluidized bed thermogravimetric analysis for chemical looping combustion and chemical looping with oxygen uncoupling. Energy fuels. 33 (1), 449–459. doi:10.1021/acs.energyfuels.8b02909

CrossRef Full Text | Google Scholar

Li, Z., Cai, J., and Liu, L. (2021). A first-principles microkinetic rate equation theory for heterogeneous reactions: application to reduction of Fe2O3 in chemical looping. Ind. Eng. Chem. Res. 60 (43), 15514–15524. doi:10.1021/acs.iecr.1c03214

CrossRef Full Text | Google Scholar

Li, Z., Zhang, W., Chen, Z., Ren, Z., Ning, S., and Li, M. (2023). A kinetics mechanism of NOX formation and reduction based on density functional theory. Sci. Total Environ. 867, 161519. doi:10.1016/j.scitotenv.2023.161519

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, L., Wang, H., Cai, J., Cai, N., and Li, Z. (2022). First-principles-based multiscale modelling of heterogeneous CoO oxidation kinetics in high-temperature thermochemical energy storage. Fuel Process. Technol. 228, 107164. doi:10.1016/j.fuproc.2022.107164

CrossRef Full Text | Google Scholar

Liu, L., Li, Z., Li, Z., Larring, Y., Li, Y., and Cai, N. (2021). Fast redox kinetics of a perovskite oxygen carrier measured using micro-fluidized bed thermogravimetric analysis. Proc. Combust. Inst. 38, 5259–5269. doi:10.1016/j.proci.2020.06.160

CrossRef Full Text | Google Scholar

Liu, T., Yang, R., Zhang, G., Wu, W., Yang, Z., Lin, R., et al. (2021). Mechanism of selective catalytic reduction of NOx with NH3 over CeO2-TiO2: insight from in-situ DRIFTS and DFT calculations. Appl. Surf. Sci. 568, 150764. doi:10.1016/j.apsusc.2021.150764

CrossRef Full Text | Google Scholar

Lu, L. (2022). GPU accelerated MFiX-DEM simulations of granular and multiphase flows. Particuology 62, 14–24. doi:10.1016/j.partic.2021.08.001

CrossRef Full Text | Google Scholar

Lu, L., Gao, X., Dietiker, J. F., Shahnam, M., and Rogers, W. A. (2021). Machine learning accelerated discrete element modeling of granular flows. Chem. Eng. Sci. 245, 116832. doi:10.1016/j.ces.2021.116832

CrossRef Full Text | Google Scholar

Lupiáñez, C., Guedea, I., Bolea, I., Díez, L. I., and Romeo, L. M. (2013). Experimental study of SO2 and NOx emissions in fluidized bed oxy-fuel combustion. Fuel Process. Technol. 106, 587–594. doi:10.1016/j.fuproc.2012.09.030

CrossRef Full Text | Google Scholar

Lyngfelt, A., Leckner, B., and Mattisson, T. (2001). A fluidized-bed combustion process with inherent CO2 separation; application of chemical-looping combustion. Chem. Eng. Sci. 56 (10), 3101–3113. doi:10.1016/S0009-2509(01)00007-0

CrossRef Full Text | Google Scholar

Madon, R. J., Braden, D., Kandoi, S., Nagel, P., Mavrikakis, M., and Dumesic, J. A. (2011). Microkinetic analysis and mechanism of the water gas shift reaction over copper catalysts. J. Catal. 281 (1), 1–11. doi:10.1016/j.jcat.2011.03.008

CrossRef Full Text | Google Scholar

Maestri, M., Vlachos, D. G., Beretta, A., Groppi, G., and Tronconi, E. (2008). Steam and dry reforming of methane on Rh: microkinetic analysis and hierarchy of kinetic models. J. Catal. 259 (2), 211–222. doi:10.1016/j.jcat.2008.08.008

CrossRef Full Text | Google Scholar

Motagamwala, A. H., and Dumesic, J. A. (2021). Microkinetic modeling: a tool for rational catalyst design. Chem. Rev. 121 (2), 1049–1076. doi:10.1021/acs.chemrev.0c00394

PubMed Abstract | CrossRef Full Text | Google Scholar

Ning, H., Tang, R., Li, C., Gu, X., Gong, Z., Zhu, C., et al. (2025). Recent advances in process and materials for dry desulfurization of industrial flue gas: an overview. Sep. Purif. Technol. 353, 128425. doi:10.1016/j.seppur.2024.128425

CrossRef Full Text | Google Scholar

Noor, T., Qi, Y., and Chen, D. (2020). Hydrogen dependence of the reaction mechanism and kinetics of water gas shift reaction on Ni catalyst: experimental and DFT study. App. Catal. B Environ. 264, 118430. doi:10.1016/j.apcatb.2019.118430

CrossRef Full Text | Google Scholar

Nørskov, J. K., Studt, F., Abild-Pedersen, F., and Bligaard, T. (2014). Fundamental concepts in heterogeneous catalysis. Hoboken, New Jersey: John Wiley and Sons, Inc.

Google Scholar

Osman, A. I., Hefny, M., Abdel Maksoud, M. I. A., Elgarahy, A. M., and Rooney, D. W. (2021). Recent advances in carbon capture storage and utilisation technologies: a review. Environ. Chem. Lett. 19, 797–849. doi:10.1007/s10311-020-01133-3

CrossRef Full Text | Google Scholar

Petersen, E. E. (1957). Reaction of porous solids. AIChE J. 3 (4), 443–448. doi:10.1002/aic.690030405

CrossRef Full Text | Google Scholar

Rawal, T. B., Le, D., Hooshmand, Z., and Rahman, T. S. (2021). Toward alcohol synthesis from CO hydrogenation on Cu(111)-supported MoS2–predictions from DFT+KMC. J. Chem. Phys. 154, 174701. doi:10.1063/5.0047835

PubMed Abstract | CrossRef Full Text | Google Scholar

Razdan, N. K., and Bhan, A. (2021). Catalytic site ensembles: a context to reexamine the Langmuir-Hinshelwood kinetic description. J. Catal. 404, 726–744. doi:10.1016/j.jcat.2021.09.016

CrossRef Full Text | Google Scholar

Reuter, K. (2016). Ab initio thermodynamics and first-principles microkinetics for surface catalysis. Catal. Lett. 146 (3), 541–563. doi:10.1007/s10562-015-1684-3

CrossRef Full Text | Google Scholar

Sakai, M., Abe, M., Shigeto, Y., Mizutani, S., Takahashi, H., Viré, A., et al. (2014). Verification and validation of a coarse grain model of the DEM in a bubbling fluidized bed. Chem. Eng. J. 244, 33–43. doi:10.1016/j.cej.2014.01.029

CrossRef Full Text | Google Scholar

Sandmann Jr, C. W., and Zygourakis, K. (1986). Evolution of pore structure during gas-solid reactions: discrete models. Chem. Eng. Sci. 41 (4), 733–739. doi:10.1016/0009-2509(86)87152-4

CrossRef Full Text | Google Scholar

Sedghkerdar, M. H., and Mahinpey, N. (2015). A modified grain model in studying the CO2 capture process with a calcium-based sorbent: a semianalytical approach. Ind. Eng. Chem. Res. 54 (3), 869–877. doi:10.1021/ie503989n

CrossRef Full Text | Google Scholar

Srinivasakannan, C., Al Shoibi, A., and Balasubramanian, N. (2012). Combined resistance bubbling bed model for drying of solids in fluidized beds. Heat. Mass Transf. 48, 621–625. doi:10.1007/s00231-011-0909-2

CrossRef Full Text | Google Scholar

Stearn, A. E., and Eyring, H. (1940). Absolute rates of solid reactions: diffusion. J. Phys. Chem. 44 (8), 955–980. doi:10.1021/j150404a001

CrossRef Full Text | Google Scholar

Szekely, J., and Evans, J. W. (1970). A structural model for gas-solid reactions with a moving boundary. Chem. Eng. Sci. 25 (6), 1091–1107. doi:10.1016/0009-2509(70)85053-9

CrossRef Full Text | Google Scholar

Szekely, J., and Evans, J. W. (1971). A structural model for gas-solid reactions with a moving boundary-II. Chem. Eng. Sci. 26 (11), 1901–1913. doi:10.1016/0009-2509(71)86033-5

CrossRef Full Text | Google Scholar

Thybaut, J. W., Sun, J., Olivier, L., Van Veen, A. C., Mirodatos, C., and Marin, G. B. (2011). Catalyst design based on microkinetic models: oxidative coupling of methane. Catal. Today 159 (1), 29–36. doi:10.1016/j.cattod.2010.09.002

CrossRef Full Text | Google Scholar

Van Belleghem, J., Yang, J., Janssens, P., Poissonnier, J., Chen, D., Marin, G. B., et al. (2022). Microkinetic model validation for Fischer-Tropsch synthesis at methanation conditions based on steady state isotopic transient kinetic analysis. J. Ind. Eng. Chem. 105, 191–209. doi:10.1016/j.jiec.2021.09.017

CrossRef Full Text | Google Scholar

Wang, H., Li, Z., Fan, X., and Cai, N. (2017). Rate-equation-based grain model for the carbonation of CaO with CO2. Energy fuels. 31 (12), 14018–14032. doi:10.1021/acs.energyfuels.7b02751

CrossRef Full Text | Google Scholar

Wang, H., Li, Z., Li, Y., and Cai, N. (2021). Reduced-order model for CaO carbonation kinetics measured using micro-fluidized bed thermogravimetric analysis. Chem. Eng. Sci. 229, 116039. doi:10.1016/j.ces.2020.116039

CrossRef Full Text | Google Scholar

Wang, S., Hu, C., Luo, K., Yu, J., and Fan, J. (2023). Multi-scale numerical simulation of fluidized beds: model applicability assessment. Particuology 80, 11–41. doi:10.1016/j.partic.2022.11.011

CrossRef Full Text | Google Scholar

Wang, R., Li, Z., and Liu, L. (2025). Multiscale modeling for the reduction kinetics of a perovskite oxygen carrier based on quantum chemistry and CFD–DEM. Carbon Capture Sci. Technol. 14, 100357. doi:10.1016/j.ccst.2024.100357

CrossRef Full Text | Google Scholar

Wang, X., Wang, W., Xiong, W., Jiang, X., Ouyang, T., Bai, Y., et al. (2024). A DFT study of the mechanism of NH3-SCR NOx reduction over Mn-doped and Mn–Ti co-doped CoAl2O4 catalysts. J. Mat. Chem. C 12, 5073–5082. doi:10.1039/D4TC00239C

CrossRef Full Text | Google Scholar

Wang, Y., Liu, H., Duan, Q., and Li, Z. (2024). Understanding the negative Apparent activation energy for Cu2O and CoO oxidation kinetics at high temperature near equilibrium. Catalysts 14 (11), 832. doi:10.3390/catal14110832

CrossRef Full Text | Google Scholar

Willis, M. D., and Wilson, K. R. (2022). Coupled interfacial and bulk kinetics govern the timescales of multiphase ozonolysis reactions. J. Phys. Chem. A 126 (30), 4991–5010. doi:10.1021/acs.jpca.2c03059

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, X., Li, Y., Zhao, J., and Wang, Z. (2020). Density functional theory study on CO2 adsorption by Ce-promoted CaO in the presence of steam. Energy fuels. 34 (5), 6197–6208. doi:10.1021/acs.energyfuels.0c00972

CrossRef Full Text | Google Scholar

Yang, W., Cloete, S., Morud, J., and Amini, S. (2016). An effective reaction rate model for gas-solid reactions with high intra-particle diffusion resistance. Int. J. Chem. React. Eng. 14 (1), 331–342. doi:10.1515/ijcre-2015-0127

CrossRef Full Text | Google Scholar

Yin, P., Yang, Y. S., Chen, L. F., Xu, M., Chen, C. Y., Zhao, X. J., et al. (2020). DFT study on the mechanism of the water gas shift reaction over NixPy catalysts: the role of P. J. Phys. Chem. C 124 (12), 6598–6610. doi:10.1021/acs.jpcc.9b11142

CrossRef Full Text | Google Scholar

Yu, Y., Ji, P., Zhang, W., Yang, K., and Zhang, M. (2023). Theoretical study coupling DFT calculations and kMC simulation of CO methanation on Ni(111) and Ni3Fe(111). New J. Chem. 47, 17923–17936. doi:10.1039/D3NJ03516F

CrossRef Full Text | Google Scholar

Glossary

Ci Concentration of lattice species i, mol/m3

Cmax Maximum concentration of lattice oxygen, mol/m3

d* Number of neighboring sites of a Mn site

d# Number of neighboring sites of a Ca site

dp Particle diameter, m

Deff Effective intraparticle diffusivity, m2/s

Dj Diffusivity of gas species j, m2/s

fdrag Drag force, N

fN Normal collision force, N

fcoll Collision force, N

fT Tangential collision force, N

g Gravitational acceleration, m/s2

h Planck constant, J·s

J Flux of oxygen on the grain surface, mol/(m2·s)

kB Boltzmann constant, J/K

k+ Rate constant of an elementary reaction (forward)

k Rate constant of an elementary reaction (reversed)

Li Li number

m Mass of oxygen carrier in the reactor, kg

mox Mass of fully-oxidized oxygen carrier, kg

mp Mass of a particle, kg

mp,ox Mass of a fully-oxidized particle, kg

mre Mass of fully-reduced oxygen carrier, kg

m˙p,j Mass transfer rate of gas species j, kg/s

MO Molar mass of oxygen atoms, kg/mol

N*,tot Number of Mn sites on a surface

N#,tot Number of Mn sites on a surface

N*#,tot Number of Mn–Ca site pairs on a surface

Ni Number of surface species i

N˙ Net reaction frequency on the surface, s−1

p Gas partial pressure field, Pa

pi Partial pressure of gas species i, Pa

ps Gas partial pressure at particle surface, Pa

qi0 Partition function of species i based on the reference energy

qTS0, Partition function of transition state excluding dissociation

r Radius coordinate, m

rg Grain radius, m

rp Particle radius, m

R Ideal gas constant, J/(mol·K)

ROC Oxygen capacity of oxygen carrier

t Time, s

T Reaction temperature, K

uf Fluid velocity, m/s

vp Particle velocity, m/s

ΔvN Normal relative particle velocity, m/s

ΔvT Tangential relative particle velocity, m/s

V Gas volume, m3

Vcell Cell volume, m3

xp Particle position, m

ΔxN Normal particle overlap, m

ΔxT Tangential particle overlap, m

X Conversion of oxygen carrier in the reactor

Xg Conversion of a grain

Xp Conversion of a particle

α Particle porosity

αf Volume fraction of fluid phase

βp Drag force coefficient, kg/(m3·s)

γN Normal damping coefficient, N·(m/s)−1

γT Tangential damping coefficient, N·(m/s)−1

δ Atomic number in CaMn0.375Ti0.5Fe0.125O3−δ

Δε0,+ Energy barrier of an elementary reaction (forward), J

Δε0, Energy barrier of an elementary reaction (reversed), J

η Effectiveness factor of intraparticle gas diffusion

θi Coverage of Mn-site surface species i

θ˙ Net reaction rate per Mn site, s−1

ϑi Coverage of Ca-site surface species i

κN Normal spring coefficient, N/m

κT Tangential spring coefficient, N/m

Λ Areal density of Mn sites, mol/m2

νeff Effective kinematic viscosity, m2/s

ρf Fluid density, kg/m3

ρj Density of gas species j, kg/m3

ρt Particle true density, kg/m3

τ Viscous stress tensor, m2/s2

ϕ Thiele modulus

ω Intrinsic rate constant for gas, s−1

Keywords: multiscale model, first-principles calculation, dual-site microkinetics, bulk diffusion, CFD-DEM, oxygen carrier, chemical looping

Citation: Wang R, Li Z and Liu L (2025) Multiscale reaction model coupling dual-site microkinetics with bulk diffusion and CFD–DEM for a perovskite oxygen carrier. Front. Chem. 13:1656180. doi: 10.3389/fchem.2025.1656180

Received: 29 June 2025; Accepted: 14 August 2025;
Published: 04 September 2025.

Edited by:

Kun Zhao, Chinese Academy of Sciences (CAS), China

Reviewed by:

Jin-Woo Kim, Seoul National University, Republic of Korea
Francesco Orsini, Polytechnic University of Turin, Italy

Copyright © 2025 Wang, Li and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Zhenshan Li, bGl6c0B0c2luZ2h1YS5lZHUuY24=

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.