Numerical modelling of ceria undergoing reduction in a particle–gas counter-flow: Effects of chemical kinetics under isothermal conditions

Journal Pre-proofs Numerical modelling of ceria undergoing reduction in a particle–gas counterflow: effects of chemical kinetics under isothermal conditions Sha Li, Vincent M. Wheeler, Apurv Kumar, Wojciech Lipiński PII: DOI: Reference:

S0009-2509(20)30085-3 https://doi.org/10.1016/j.ces.2020.115553 CES 115553

To appear in:

Chemical Engineering Science

Received Date: Revised Date: Accepted Date:

6 November 2019 19 January 2020 5 February 2020

Please cite this article as: S. Li, V.M. Wheeler, A. Kumar, W. Lipiński, Numerical modelling of ceria undergoing reduction in a particle–gas counter-flow: effects of chemical kinetics under isothermal conditions, Chemical Engineering Science (2020), doi: https://doi.org/10.1016/j.ces.2020.115553

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier Ltd.

Numerical modelling of ceria undergoing reduction in a particle–gas counter-flow: effects of chemical kinetics under isothermal conditions Sha Lia, Vincent M. Wheelerb, Apurv Kumara, c, and Wojciech Lipińskia,* a

Research School of Electrical, Energy and Materials Engineering, The Australian National University, Canberra,

ACT 2601, Australia b

Department of Engineering and Technology, University of Wisconsin–Stout, Menomonie, WI 54751, USA

c

School of Science, Engineering and Information Technology, Federation University, Mt. Helen, Ballarat, Australia

Highlights 

Particle–gas reactive counter-flow analyzed using Euler–Lagrange approach;



Factors limiting the outlet reduction extent identified at varying conditions;



Design strategies proposed to approach the upper limit of reaction extent;



Correlation developed relating reduction extent to flow rate and kinetics.

Abstract A numerical model is employed to simulate a single tube reactor featuring a downward particle flow counter to an upward inert gas flow for ceria reduction in the dilute flow regime. The coupled phenomena of mass and momentum transfer as well as chemical kinetics are simulated assuming isothermal operation for the reactor. The model predicts the reduction extent under varying reaction kinetics as well as design and operational choices. The reduction extent is found to increase with the reaction rate constant until achieving the thermodynamic upper limit at a certain critical value. This critical rate constant signifies a transition from a chemical kinetics limited conversion to a gas advection limited conversion. The effect of the reactor length and the particle size on reaction extent is studied for a range of realistic cases. An empirical correlation is

developed to quantify the effects of particle and gas flow rates on reduction extent at both slow and fast kinetics. The present work offers insights to help guide reactor design and operation towards achieving the maximum reduction extent. Keywords: particle–gas reactive counter-flow; Euler–Lagrange approach; ceria reduction; reaction kinetics; thermodynamics 1. Introduction Solar-driven, non-stoichiometric, ceria-based reduction/oxidation (redox) cycling for splitting of water and/or carbon dioxide is a promising thermochemical pathway for synthetic hydrocarbon fuel production [1, 2]. The ceria-based redox cycle for water splitting is composed of an endothermic reduction step, 1 1 1 CeO 2 ox (s)  CeO 2 red (s)  O 2 (g),   2

(1)

and an exothermic oxidation step, 1 1 CeO 2 red (s)  H 2O(g)  CeO 2 ox (s)  H 2 (g).  

(2)

Fuel production is proportional to the non-stoichiometry change of ceria which can be driven by temperature-swing [3] or pressure-swing [4, 5] of the reduction/oxidation chamber(s). A variety of reactor concepts/design have been analyzed/tested to implement ceria-based redox cycling, including: semi-batch operation using a fixed bed [1, 2]; continuous operation using a solid–gas co-current/parallel flow (PF) or countercurrent flow (CF) configurations [6–8]; and continuous operation using an isothermal membrane [5, 9], among others. Our previous thermodynamic analyses [10, 11] have identified the upper limit of reduction extent in a generic reactor implementing a solid–gas CF configuration with prescribed inlet flow rates and

thermodynamic states. However, the practical realization of the CF configuration is difficult. The design of a solid–gas CF reactor critically relies on the implementation of moving solids using either rotary structures [8] or a downward flow of particles falling by gravity [6, 7]. Particle flow reactors (or aerosol flow reactors) have found wide application in many solar thermochemical processes [6, 7, 12–14] due to the following advantages as compared to fixed-bed reactors [6, 15]: (i) reactor operation is continuous, and (ii) heat and mass transfer rates are high due to high specific surface area at the particle–gas interface. A typical particle flow reactor usually consists of multiple tubes containing the particle–gas flow irradiated by concentrated solar energy to allow for high particle throughput and efficient solar utilization [16–18]. For a given tubular reactor where particles are subject to a short residence time, the reduction extent relies strongly on the coupled phenomena of mass, momentum and heat transfer as well as chemical kinetics. The early reactive flow models in solar thermochemistry are mainly concerned with radiation–chemistry interactions, while hydrodynamics is trivial or not considered at all [15, 19–21]; see the review article by Lipiński et al. [22]. For those that do consider the hydrodynamic characteristics, most are focused on the downward PF configuration [16, 23, 24] and the CF configuration has been largely neglected [7]. A widely employed assumption for the PF configuration is perfect entrainment based on low Stokes numbers (<<1) when using very small particles [14, 16, 23–25], though treatment of velocity non-equilibrium can also be found [26, 27]. Unlike the PF configuration, the CF configuration no longer warrants the particle entrainment assumption and hence the particle–gas hydrodynamic interaction must be resolved using separate momentum equations for each phase. This can be done using either the Euler–Lagrange approach to track all particles (or their representative set) or the Euler–Euler approach to treat both phases as interpenetrating continua [28]. Typically, the former approach is more applicable for dilute

particle flow regime while the latter for dense flow regime, though application of the Euler–Euler approach to dilute flow can also be found [16, 23, 24]. Ceria reduction using a CF reactor has earlier been studied experimentally and numerically. Welte et al. demonstrated the performance of an aerosol reactor prototype implementing a CF configuration under realistic concentrated radiation [7] and controlled heating rate [6] for the first time. In addition, a coupled heat and mass transfer model was also formulated to help elucidate the conversion-limiting mechanisms [7]. They concluded that the reaction extent is mainly limited by gas advection for smaller particles and by heat transfer for larger particles and/or higher particle flow rates. However, reaction kinetics is not considered as a limiting factor due to the assumption of chemical equilibrium everywhere in their model [7]. Though ceria is widely recognized to offer fast kinetics [1], achieving a certain equilibrium state still needs a finite time, which may be similar to or longer than the particle residence time of a flow reactor—typically several seconds [6, 23]. Consequently, it is possible for the reaction kinetics to become the conversion-limiting factor and therefore a realistic kinetic model should be employed to test this possibility. Furthermore, heat transfer rate plays an important role in determining the reaction extent and may become the conversion-limiting factor too, particularly if large temperature gradients exist within a reactor [1]. Previous numerical studies on flow reactors have indicated that both phases are rapidly heated from ambient condition to the same maximum temperature that can be maintained uniformly over a long reactor zone when particle flow is dilute and particle size is small [7, 23, 24]; heat transfer is unlikely to limit the reduction extent under such conditions [7]. Our previous thermodynamic analyses [10, 11] were limited to an abstract representation of a thermochemical reactor where a reduction extent upper limit was determined for ceria for given inlet flow rates and thermodynamics states. Mass and momentum transfer as well as reaction

kinetics were not considered [10, 11]. Herein, we offer an enhanced understanding of the particle– gas CF configuration for realistic systems. A single tube reactor is adopted as the model system composed of a downward particle flow against an upward gas flow at isothermal conditions. The ceria reduction step is studied in detail and the oxidation step is left for further work. A Euler– Lagrange approach is employed to model the reactor in the dilute particle flow regime. The coupled phenomena of mass and momentum transfer as well as chemical kinetics are captured. Particles with uniform size below 100 μm are selected in our numerical simulation to enable the simplifications described in Section 2. The effects of selected operational and design parameters will be investigated with the following objectives: (i) to identify the conversion-limiting factor at varying kinetics and design choices; (ii) to understand the effects of operating conditions on reaction extent and to see how much it deviates from its thermodynamic limit; (iii) to approach this reduction upper limit via proper reactor design and operation. 2. Problem statement A schematic of the model system is illustrated in Fig. 1(a) for a two-phase domain composed of a downward ceria particle flow counter to an upward gas flow within a single vertical tube reduction reactor. Inert nitrogen flow is supplied to the reactor from the gas inlet at the bottom with a mass flow rate of mg,in and an oxygen partial pressure of pO2 ,in . Ceria particles are fed from the opposite end with a mass flow rate of mp,in distributed uniformly over the inlet cross-section area at an oxidized state of  ox .

(a)

(b)

Fig. 1. Schematic of (a) the model system assuming isothermal operation; (b) four steps involved in the interphase oxygen transport process for a reacting spherical particle: (i) lattice diffusion within the particle volume, (ii) surface reaction at the particle–gas interface using doubly ionized oxygen vacancy mechanism as an example, (iii) oxygen diffusion within the gas film, and (iv) bulk gas advection of product oxygen towards the main gas flow.

The coupled phenomena under consideration for the isothermal multiphase flow system include momentum and mass transfer as well as chemical kinetics. The injected particles must be subjected to a higher graviational pull than the upward forces from the counter flowing gas to realize the particle–gas CF configuration. Thus, the gas flow rate and ceria particle diameter dictate possibility of achieving a CF configuration in the reactor. During particle-gas contact, the particles undergo thermal reduction and release oxygen at the particle surface that can be swept away by the gas flow, causing interphase mass transfer. There are four steps involved in the interphase oxygen transport process for a reacting ceria particle as illustrated in Fig. 1(b): (i) lattice diffusion of solid oxygen within the particle volume towards the particle surface; (ii) surface reaction at the particle–gas interface; (iii) film diffusion of product oxygen from the particle surface towards the

bulk gas; and (iv) bulk gas advection of product oxygen towards the main gas flow. This transferred oxygen will serve as a mass source to the gas flow, leading to an increase in both oxygen concentration and gas flow rate along the flow path. On the other hand, the non-stoichiometry state of the particle phase will become more reduced as the reaction proceeds. Consequently, a steadystate distribution of oxygen partial pressure and particle non-stoichiometry state will be established along the flow path. The gas flow exits the reactor with an oxygen partial pressure of pO2 ,out while the particle flow leaves at an reduced non-stoichiometry state of  red . A number of assumptions are made to assist in the model development: (i) the model system is operated at steady-state assuming isothermal condition: Tred = Tp = Tg ; (ii) the gas flow is assumed to be laminar and fully developed, and the particle phase is in the dilute flow regime; justification on this assumption can be found in section 4.3; (iii) the effects of particle sintering and aggolomeration are not considered; (iv) all particles are assumed non-porous, spherical, and uniform in size, i.e., the effect of particle size distribution is not considered; (v) the particle diameter is assumed to remain unchanged during the reduction process, so the particle density will decrease accordingly; (vi) lattice- and film-diffusion steps are assumed fast enough as compared to the suface reaction step such that each single particle has uniform composition throughout and the gas concentration at the particle surface is identical to that in the bulk gas; justification of this assumption is detailed in the Supplementary Information for interested readers. (vii) all particles are rigid and non-rotating, so each particle has a single velocity; (viii) steady-state drag is the only significant fluid dynamic force acting on the particles; (ix) the efflux velocity of product oxygen is assumed uniform over the particle surface and the blowing effect of product oxygen from particle surface on the drag force is not considered [28]; (x) the particle volume fraction is so small such that the gas volume fraction is assumed to be unity in gas phase equations.

3. Two-phase reactive flow model A two-way coupled Euler–Lagrange approach is employed to model the dilute particle–gas CF configuration in the tube. The gas phase is treated as a continuum by solving Eulerian volumeaveraged governing equations, while the discrete particle phase is studied by tracking all particles in the flow field within the reactor. 3.1 Particle-phase equations The underlying physics of the discrete particle phase is described by ordinary differential equations (ODEs) based on a Lagrangian-tracking approach—the particle flow is resolved down to the particle level. The mass change of a reacting particle mp with constant size due to surface reaction reads dmp dt

 Vp

d p dt

  Ap M O2 rO2  Vp M O2 rO2 ,

(3)

where M represents a molar mass and rO2 as well as rO2 are the molar oxygen production rate for ceria reduction, Eq. (1), per particle surface area Ap and particle volume Vp , respectively:

rO2 

1 dnO2 1 dnO2 , and rO2  . Ap dt Vp dt

(4)

The negative sign in Eq. (3) indicates a mass sink. The particle density  p corresponds to the nonstoichiometry state  of a particle via 

p  1   

M O2    CeO2 , 2 M CeO2 

where CeO2 is the density of stoichiometric ceria. Combining Eqs. (3) and (5) yields

(5)

2M CeO2 d 12M CeO2  rO2  r  , dt d p CeO2 CeO2 O2

(6)

where d p is the particle diameter. Several existing kinetic models for ceria reduction in the form of rO2 or rO2 will be elaborated later in Section 3.3. The equations of motion for a single particle are: dxp,i dt

mp

dvp,i dt

 vp,i ,

(7)

g

 mp gi  g Vp gi  mp f drag

v

g

g,i

p

 vp,i

,

(8)

where  g is the gas density, vp,i and vg,i are the particle and gas velocities in the i direction, respectively. The angle brackets with a superscript g represent intrinsic volume-averaged properties over the gas phase; this will be elaborated later in Section 3.2. The respective force terms on the right-hand side of Eq. (8) are gravity force, buoyant force, and steady-state drag force. Particle–particle interactions are not considered here since these effects are negligible in the dilute flow regime. The symbols  p and f drag appearing in Eq. (8) are the particle relaxation time (or particle response time) and the drag factor, respectively, as expressed by:

p dp2 p  , 18g f drag 

(9)

Cdrag Rer 24

,

(10)

where  g is the gas phase viscosity and Re r is the relative Reynolds number as defined by:

Rer 

g

g

d p vg,i

g

g

 vp,i .

(11)

The drag coefficient Cdrag appearing in Eq. (10) can be determined by several correlations available in literature as a function of the relative Reynolds number Re r [28–32]. The spherical law by Morsi and Alexander [32] is adopted here and will not be detailed. Integration of the coupled Equations (3), (7) and (8) over time for all particles will yield the particle mass (also density) and velocity in the solution domain. Unlike the gas flow that has a continuous distribution in the field, the particle flow is discretized into a number of particle streams that satisfy mp 

N total

m k 1

p,k

,

(12)

where mp and mp,k are the total particle mass flow rate and the mass flow rate of particle stream k, respectively, and N total is the total number of particle streams. In a steady flow calculation, each particle represents a stream of many particles that follow the same trajectory. The relationship between mass flow rate mp,k and particle mass mp,k associated with a particle stream k is given by

mp,k  Np,k mp,k ,

(13)

where Np,k is the number flow rate of particle stream k. The properties associated with each single particle ( mp ,  p and vp,i ) as well as with each particle stream k ( mp,k and Np,k ) therefore offer a complete description of the particle flow. This provides a basis for determining the accumulated discrete particle model (DPM) source terms and particle volume fraction within a computational cell in the flow field used for developing the volume-averaged gas phase equations. The accumulated particle mass in a given cell Vc due to particle surface reaction is expressed per volume of the cell by examining the change in mass flow rate of all particle streams traversing the cell:

S DPM,mass 

1 Vc

 (m

p,in,k

 mp,out,k ) 

k

1 Vc

N

p,k

(mp,in,k  mp,out,k )

(14)

k

The mass exchange term S DPM,mass will appear as a source term in the gas phase continuity, Eq. (18), and oxygen species conservation, Eq. (19). Similarly, the accumulated particle momentum within a computational cell in the field is expressed as:

SDPM,momentum,i

1  Vc

g     g   k Np,k  mp,in,k vp,in,k ,i  mp,out,k vp,out,k ,i   mp,out,k   gi   gi  tk  , (15) p    

where tk is the residence time of particle stream k within the cell and the summation is carried out for all particle streams traversing the cell. The subscript i represents Cartesian coordinate direction. Note that the external forces of particle gravity force and buoyancy force do not contribute to the momentum source of the gas phase and must be excluded; see the second term in the square bracket of Eq. (15). This momentum exchange term S DPM,momentum,i will appear as a source term in the gas phase momentum equation (see Eq. (21)) as will be described in section 3.2. The local particle volume fraction within a computational cell Vc is calculated by

f v,p 

Vp N p Vc



Vp  N p,k tk k

Vc

.

(16)

where N p is the number of particles in a computational cell Vc . So, the local gas phase volume fraction reads f v,g  1  f v,p ,

(17)

which will also appear in the gas phase governing equations arising from volume-averaging process in the coming section. 3.2 Gas-phase equations

The gas phase equations governing mass and momentum transfer as well as chemical kinetics are developed based on the volume-averaging approach [33]. Unlike the Lagrangian particle phase equations that are resolved down to the particle level, the Eulerian gas phase equations are based on discretizing the flow field into computational cells that are essentially averaging volumes containing a group of particles. This allows the information about transport phenomena at the particle level to be communicated to the macroscale by introducing effective transport properties that capture the influence of the particles. The steady-state mass conservation equation for the gas phase reads



g

 f v,g g

g

vg,i

xi

S

DPM,mass

.

(18)

Please note that the Einstein convention is used here for all gas phase equations where double appearances of indices imply summation. The steady-state oxygen species conservation accounts for bulk advection and species diffusion as well as the interphase mass transfer from particles due to surface reduction:



 f v,g YO2

g

g

xi

where YO 2

g

g

vg,i

g

    D xi  

O2 ,eff



 YO2

g

g

xi

g

   S   

DPM,mass

,

(19)

is the volume-averaged mass fraction for oxygen species in the gas phase. The

effective diffusion coefficient DO2 ,eff is determined following a similar logic to determine the effective thermal conductivity in the work by Crowe et al. [28]: DO2 ,eff  f v,g DO2 ,g  f v,p DO,p ,

(20)

where DO2,g is the mass diffusion coefficient of oxygen in the gas phase [34] and DO,p is the ambipolar diffusion coefficient of oxygen within the particle [35].

The steady-state momentum conservation for the gas phase is given by:



 f v,g g

g

vg,i

g

vg, j

g

x j   v   g ,i  f v,g g   x j x j   

g



 f

g

v,g

 vg , j xi

g

g

gi  f v,g

 p xi

g

  vg , m   f v,g  g,bulk  2 g     3  xm  

g

  ij   SDPM,momentum,i  

(21)

where  g and  g,bulk are the dynamic viscosity and the coefficient of bulk viscosity for the gas phase, and  ij is the isotropic tensor Kronecker delta. Note that Stokes’ hypothesis is adopted here to treat the coefficient of bulk viscosity  g,bulk to be zero [36]. The viscosity of gas mixture  g is computed using a simple mole-fraction-weighted mixing rule by Graham [37]:  Y g i g    i Mi i  O2 ,N 2  



Yj

j  O2 ,N 2

g

Mj

 ,  

(22)

This mixing rule is selected because it is adequate for gas mixtures in which the components have nearly the same molecular weight [38], like the components of nitrogen and oxygen being studied here. The species viscosities are evaluated as functions of temperature [39]. The ideal gas law is employed to offer a closure relationship between pressure and gas phase density,

p

g

g  Yi       i O2 ,N2 M i  g  

g

g

R Tg .

(23)

Note that the gas phase volume fraction appearing in Eqs. (18) (19) and (21) is approximated as 1 based on assumption (x) in Section 2. 3.3 Reaction kinetics Several models have been proposed for the rate expressions, rO2 or rO2 , appearing in Eq. (3) for ceria reduction. Oles and Jackson [27] developed a comprehensive kinetic model considering the

Table 1. Summary of reaction rate models by Keene et al. [40, 41] and Bulfin et al. [42] Reaction mechanism

CeCe  OO

 kred,K 2  kox,K 2

2CeCe  OO

 kred,K 3

2CeCe  OO

 kred,Bulfin

 kox,K 3

 kox,Bulfin

1 CeCe  VO  O2 (g) 2 1 2CeCe  VO  O2 (g) 2

1 2CeCe  VO  O2 (g) 2

General rate expression

Rate

Ref(s)

 2 CCe CO  kox,K  2 CCeCe CV pO1/22 rO2  kred,K

(24)

constants unknown

[41]

2 2 1/2  3 CCe  3 CCe rO2  kred,K  C   kox,K Ce CV pO2 O

(25)

unknown

[15, 40]

  rO2  kred,Bulfin CO  kox,Bulfin CV pOn 2

(26)

known

[42]

Ce

O

Ce

O

O

O

O

O

lattice diffusion step, the surface exchange step and the oxygen desorption step on the particle level to account for the concentration gradient within each single particle. Such detailed information is certainly valuable in predicting the reactor performance. However, as a first-stage analysis, our focus is on the overall particle phase behavior, and a simplification is thus made to treat each particle with uniform composition throughout (see assumption (vi) in Section 2), so their model is not adopted here. Keene et al. [40, 41] developed two mechanistic models assuming either a singlyor doubly-ionized oxygen vacancy defect mechanism based on the law of mass action, crystal site conservation and electroneutrality. However, due to a lack of experimental data at non-equilibrium conditions, the rate constant in their models remains unknown, so equilibrium chemistry was assumed by using an extremely large rate constant [41]. Later, Bader et al. [15] adopted their model [40] to study heat and mass transfer phenomena in a particle suspension considering a wide range of rate constants to understand their effects. Another kinetic model was developed by Bulfin et al. [42] in the form of a reversible Arrhenius equation assuming a constant site fraction for the removable oxygen within ceria in the regime beyond dilute species reaction. The model parameters are fixed by fitting experimental data at equilibrium from literature and at non-equilibrium from their fixed-bed measurement. Table 1 summarizes the models developed by Keene et al. [40, 41]

and by Bulfin et al. [42] using general rate expressions to better showcase their fundamental differences regarding the reaction mechanism. The detailed derivation of Eqs. (24)–(26) will not be elaborated here and readers can consult the original works [40–42]. The final versions of the two models by Keene et al. [40, 41], Eq. (24) and (25), are written here in terms of d dt using to Eq. (6):

d dt d dt

(27)

Keene,VO

1/2  2  2  hO02 ( )  Tp sO02 ( )  2  pO2   6M CeO2 nsf,K  exp        ,  CeO2 d p  2 RT p p ref      

(28)

Keene,VO

1/2  3  3  hO02 ( )  Tp sO02 ( )  3  pO2   6M CeO2 nsf,K  exp        ,  CeO2 d p  2 RTp pref       

 2 and nsf,K  3 are the undetermined rate constants for the singly- and doubly-ionized where nsf,K

oxygen vacancy mechanisms, respectively. The model by Bulfin et al. [42], Eq. (26), written in terms of the non-stoichiometry coefficient is:

d dt

Bulfin

   N Bulfin  kred,Bulfin ( x   )  kred,Bulfin  pOn 2 

(29)

where we introduce a dimensionless rate modification factor NBulfin to account for kinetics uncertainties associated with the particle flow reactor as will be elaborated in the Supplementary Information. The rate constants in Eq. (29) are expressed in the Arrhenius format,

E  ki,Bulfin  Ai exp  a,i  , i  red,ox .  RT   p

(30)

Next, we need to decide which kinetic model to use. First, the two models by Keene et al. [40, 41] were developed based on the equilibrium experimental data by Panlener et al. [43], as is consistent with our previous thermodynamic models [10, 11]. By contrast, the kinetic model by Bulfin et al. [42] predicts slightly different equilibrium conditions as confirmed by the authors

themselves in another work [44]. Second, the doubly-ionized oxygen vacancy ( VO ) mechanism is more widely accepted in the field as compared to the first reaction mechanism ( VO ) in Table 1 [35, 43, 45]. Therefore, the VO model by Keene et al. (Eq. (28)) is selected in the present work. 3.4 Initial and boundary conditions The particle phase equations are in the form of ODEs, so initial conditions are required for the particle flow to specify initial particle mass as well as velocity at injection: mp

vp

tp  0

tp  0



 6

d p3  p,in

0   0    vp,z ,in

  .  

(31)

(32)

When a particle reaches a physical boundary, a condition is also required to determine the fate of the particle trajectory at that boundary. The condition at the tube wall is set as an inelastic rebound type with a loss in momentum by defining normal ( enorm ) and tangential ( etangent ) coefficients of restitution [46], ew  max(0.7,1  0.0136 impact ) , ( w  norm, tangent)

(33)

to relate the particle normal and tangential velocities before and after collision v1p,w

x 2  y 2  Rt,in

 ew vp,0 w

x 2  y 2  Rt,in

, ( w  norm, tangent)

(34)

though such a likelihood is quite slim due to the dilute particle flow regime under consideration. The symbol appearing in Eq. (33),  impact , is the particle impact angle defined as the acute angle between particle velocity vector and wall tangential vector. The subscripts of norm and tangent indicate normal and tangential directions, respectively, and the superscripts of 0 and 1 indicate

before and after collision. All particle properties after collision will remain the same as before collision. The particle tracking will terminate when the particle hits the outlet boundary. The gas phase equations are in the form of partial differential equations (PDEs) under steady state, so only boundary conditions are required. A prescribed mass flow rate and oxygen species mass fraction is applied to the gas inlet boundary. A no-slip condition with zero diffusive species flux is specified at the tube wall. The boundary condition at the gas outlet is prescribed with an ambient pressure. Other variables at outlet can be extrapolated from solutions within the interior domain. 3.5 Numerical solution To numerically solve the coupled, multiphase and multiphysics 3D governing equations, the commercial software package ANSYS Fluent 17.1 is employed [47]. Coupling between particle and gas phases relies on updating the DPM source terms—Eqs. (14) and (15)—appearing in the gas phase governing equations—Eqs. (18), (19) and (21)—during every gas phase solution iteration. The calculations are performed iteratively between individual phase computations until the convergence criterion is met for both phases (all residuals below 10 -6 ). Given that the underlying physics of particle phase is described by ODEs while that of gas phase by PDEs, each phase has its own numerical mechanisms and discretization schemes. The time step tp used to integrate the particle ODEs is controlled by setting a step length factor ( p,cell ) roughly equivalent to the number of steps required for a particle to traverse a gas phase computational cell. This enables capturing the detailed behavior of each single particle within a gas cell in order to get a meaningful statistical average of DPM source terms; see Eqs. (14) and (15). An explicit Euler scheme is employed to integrate the particle mass equation, Eq. (3). The equations of motion are integrated using a semi-implicit trapezoidal scheme for Eq. (7)

and a fully implicit Euler scheme for Eq. (8). A segregated algorithm for the ODE solver is used to solve the coupled mass and momentum particle phase equations sequentially until all solution convergences to a residual below 106 . Spatial discretization of the cylindrical gas domain is performed by creating structured hexahedral meshes in ICEM CFD 17.1 using an O-type grid multiblock strategy as recommended here: [48, 49]. The finite volume method (FVM) is employed to integrate the gas phase PDEs over cell volumes to construct the discrete algebraic equations. To linearize the discrete equations, different interpolation schemes are required to evaluate pertinent variables at cell faces since Fluent stores solution variables at cell centers by default. The central difference scheme and second-order upwind scheme are selected to evaluate the diffusion and advection terms in all equations, respectively. Face values of the pressure term appearing in the momentum equation are predicted using a second order interpolation scheme. To solve the resultant linear algebraic equations, a point implicit linear equation solver (Gauss–Seidel) is used in conjunction with an algebraic multigrid method. A pressure-based segregated algorithm of Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) [50] is employed to solve the coupled mass and momentum equations sequentially until convergence criteria are satisfied. The solution approach has been validated by comparing with existing results in terms of hydrodynamic and chemical behaviors after a mesh independence test, and it can be found in the Supplementary Information. 4. Results and discussion A summary of all model parameters used to investigate the effects of selected operational and design parameters under varying reaction rate constants is found in Table 2. The parametric study is performed by varying one single parameter within its range with all other parameters held

constant at their baseline values unless stated otherwise. Note that the intrinsic volume-averaged non-stoichiometry results for the particle phase ( 

p

) within a given cell volume Vc will be

presented in the following subsections in order to be consistent with the gas phase results. Table 2. Model parameters for numerical analyses Parameter

Baseline value(s)

Values for parametric study

Unit

Dt,in psys

0.0191

—

m

1

—

atm

1 10 0 7220 1 vp,ter (d p )  2mg,in At-1 g,in

—

atm

— —

— kg m-3

—

m s-1

0.35 100

[0.1, 1.2] [50, 100]

m μm

2.5 106

[ 2.5 106 , 1.0  105 ]

kg s-1

2.5 105 1900

[ 2.5 105 , 1.5  104 ] [1673, 1900]

kg s-1

1.0  10

[ 5.0  10 , 1.0  10 ]

mol m-2s-1

pO2 ,in δox ρp,in vp,z,in Lt dp mg,in

mp,in Tred  3 nsf,K

5

4

3

6

K

The baseline dimensions of the isothermal cylindrical tube are 0.0191m inner diameter and 0.35 m axial length to match an experimental setup employed by Scheffe et al. [6, 7] for the same application. The range of reactor length is chosen based on experimental setups reported in literature [6, 12] to allow for both short and long high-temperature zones. The system is operated at atmospheric pressure and the inlet oxygen partial pressure is set to be 105 atm based on available air separation technologies [51]. Particles at injection are assumed to be at a fully oxidized state in order to drive ceria reduction effectively. As a result, the particle initial density is set equal to its true material density based on the non-porous particle assumption (iv) as stated in Section 2. Choices on particle size and gas mass flow rate are closely coupled in order to realize the CF configuration. Additionally, the range of particle diameter is chosen to ensure the particle

residence time below 3.5 s; its upper limit is set as 100 μm to approximate assumption (i) [7] and assumption (vi) [35] as made in Section 2. The particle mass flow rate and particle velocity at injection are chosen such that they satisfy the dilute flow assumption (ii); the upper limit of particle mass flow rate is selected to be 175 mg s-1 in order to avoid the heat transfer-limited region as suggested by Scheffe et al. [7], thus approximating the isothermal assumption. The operating temperature is selected in the range of 1673–1900 K as typically conducted in lab experiments [1, 6, 7]. The kinetic rate constant appearing in the model by Keene et al. (see Eq. (28)) is a variable with a wide range of values being considered to allow for both slow and fast reaction kinetics. In the following subsections, the effects of reaction rate constant and selected design parameters are studied first in order to identify the conversion-limiting factors at different conditions and to help guide reactor design towards achieving higher reduction extent. Next, the effects of selected operational parameters on outlet reduction extent are investigated, with a correlation further developed to quantify their influences. 4.1 Effect of reaction rate constant For the baseline case the particles are subject to a fixed residence time below 1s, which is short enough that the final reduction extent can be strongly affected by the reaction kinetics. Fig. 2 shows the longitudinal section contour plots of oxygen mass fraction (Fig. 2(a)) and particle density (Fig.2(b)) under varying reaction rate constants based on the doubly-ionized oxygen vacancy  3 values, higher oxygen mass fraction and lower model by Keene et al. [40]. For all cases of nsf,K

particle density are observed at the centerline than at the wall zone, indicating higher reduction extents at the centerline because of the slower particle velocity hence longer residence time there. Radial gradient of oxygen mass fraction is less distinct than that of particle density at lower

reaction rate constants; this is due to the gas diffusion effect which effectively homogenizes the  3 further increases, the radial gradient in radial oxygen concentration (see Eq. (19)). When nsf,K

particle density becomes less apparent because reaction rates at both the centerline and the wall are so fast that their difference is small.

(a)

(b)

Fig. 2. Longitudinal section contour plots of (a) oxygen mass fraction in the gas phase, and (b) particle  3  1  104 , 5 104 , 1  105 , 5 105 and 1 106 density at a variety of reaction rate constants of nsf,K

molm-2 s-1 , respectively (from left to right) for the baseline case.

Given the large length-to-diameter ratio of the reactor tube, profile differences of solution variables in the axial direction are much more distinct than those occurring in the radial direction and hence are our main focus. Fig. 3 displays the steady-state centerline profiles of oxygen

(b)

(b)

(c) Fig. 3. Effect of reaction rate constant on steady-state centerline profiles of (a) oxygen partial pressure, (b) particle non-stoichiometry state, (c) local reaction rate along z direction for the baseline base; the equilibrium non-stoichiometry states corresponding to the local oxygen partial pressure is illustrated in Fig. 3(b) to facilitate our discussion.

partial pressure (Fig. 3(a)), particle non-stoichiometry state (Fig. 3(b)) and normalized local reaction rate (Fig. 3(c)) for the baseline case under varying rate constants. The equilibrium nonstoichiometry states corresponding to the local oxygen partial pressures are illustrated in Fig. 3(b) to aid our analysis. Similar results can be produced using the kinetic model by Bulfin et al. [42] and will not be detailed here. Interested readers are referred to the Supplementary Information for this information. For all rate constants, both the oxygen partial pressure and particle non-

stoichiometry increase along their respective flow paths due to the CF configuration. At lower rate constants, the increase in both pO2

g

and 

p

is slow and almost linear due to the extremely low

reaction rates everywhere as illustrated in Fig. 3(c). As a result, the reduction process is proceeding far from the equilibrium condition as shown in Fig. 3(b). When the rate constant grows larger, a change in shape/slope is observed for both pO2

g

and 

p

profiles, with much greater increase

occurring at near-inlet and -outlet boundaries than that at intermediate zones. Similar numerical profiles have been reported before [7, 52] under the assumption of equilibrium everywhere. This is consistent with the higher local reaction rates at near-inlet and -outlet regions than at intermediate zones as observed in Fig. 3(c). The behavior of local reaction rate near the particle

 3  5 104 mol m-2s-1 is fundamentally due to the reaction rate inlet region (z=0.35m) when nsf,K expression scaled with  3 (see Eq. (28)): the particles at injection is assumed at a purely oxidized state (  ox =0), so the reaction rate at z=0.35m is almost zero. As soon as the particle enters the flow,

 becomes positive and the reaction rate starts to increase; the later decline in reaction rate observed downstream is a result of decreased reaction driving force at relatively high pO2



p

g

and

values since the particles are closer to their thermodynamic equilibrium state. Fig. 4 displays the effect of reaction rate constant on average reduction extent at the particle

outlet  red

p ave

as defined by Scheffe et al. [6, 7]:

 red

p ave



2nO2 np,in



2M CeO2 (mp,in  mp,out ) M O2

mp,in

(35)

where nO 2 is the change of oxygen molar flow rate from gas inlet to outlet, and np,in is the particle molar flow rate at entrance. The thermodynamic prediction is also included in Fig. 4 to facilitate

 3 ,cr seems to exist above which the average outlet our analysis. A critical reaction rate constant nsf,K

reduction extent does not continue to increase and always converges to the thermodynamic  3 ,cr , a higher rate constant means larger prediction. When the reaction rate constant is below nsf,K  3 ,cr signifies a transition from outlet reduction extent. Therefore, this critical rate constant nsf,K  3  nsf,K  3 ,cr to gas advection limited conversion when chemical kinetics limited conversion when nsf,K

 3  nsf,K  3 ,cr . nsf,K

Fig. 4. Effect of reaction rate constant on average reduction extent of the particle phase at exit for the baseline case; the thermodynamic prediction is added to facilitate the analysis.

4.2 Effect of selected design parameters Next, the numerical model is employed to investigate the effect of selected design parameters, reactor tube length ( Lt ) and the particle diameter ( d p ). Realistic reaction kinetics are not infinitely fast, so a short reaction time and/or slow reaction kinetics can become the conversion-limiting factor under certain design choices. Understanding their influence could offer insights that help

(a)

(b)

(c)

(d)

 3 = 2  104 mol m-2s-1 (left Fig. 5. Effect of reactor tube length on (a) average particle residence time at nsf,K y-axis) and average reduction extent at varying reaction rate constants (right y-axis), and centerline profiles of (b) oxygen partial pressure, (c) particle non-stoichiometry state, (d) local reaction rate along z-direction

 3 = 2  104 mol m-2s-1 for the baseline case; the thermodynamic prediction and at the rate constant of nsf,K local equilibrium conditions are added to facilitate our analysis.

 3 is guide reactor design towards achieving the maximum reduction extent if the rate constant nsf,K

relatively low. Fig.5(a) shows the effect of reactor tube length on average outlet reduction extent  red

p ave

at

 3 = 2  104 mol mvarying rate constants (right y-axis) and on average particle residence time at nsf,K

2 -1

s (left y-axis). Changing the reactor length will only affect the particle residence time—see the

 3 , there linear increase in particle residence time with larger Lt in Fig. 5(a). For each case of nsf,K

exists a specific critical reactor length Lt,cr above which the reduction extent does not continue to increase. This suggests that the reactor length by design should approach Lt,cr from the lower side in order to avoid wasting more energy in maintaining the isothermal condition for longer tubes. In addition, this critical length is found to decrease with higher reaction rate constant because particles will need less time thus shorter reactor length to achieve the thermodynamic reduction extent at faster reaction kinetics. Figs. 5(b)–(d) show the effect of reactor length on centerline profiles of oxygen partial pressure, particle non-stoichiometry state and local reaction rate,  3 = 2  104 mol m-2s-1 as a reference. With Lt increasing, all three respectively, using the case of nsf,K

profiles become more complete and simply expand to a longer reactor without changing their shape and magnitude when Lt is higher than Lt,cr . Fig. 6(a) shows the effect of particle diameter on average outlet reduction extent at varying  3 = 1 10 4 mol m-2s-1 (left y-axis). rate constants (right y axis) and on particle residence time at nsf,K

Figs. 6(b)–(d) display the corresponding centerline profiles of oxygen partial pressure, particle  3 = 1 10 4 mol m-2s-1; these non-stoichiometry state and local reaction rate, respectively, at nsf,K

profiles are quite similar to those in Fig. 3 under the effect of reaction rate constant and hence will  3 , a certain critical particle diameter d p,cr seems to exist, not be detailed here. For all cases of nsf,K

below which there is no further gain in the average outlet reduction extent. The larger reduction extent observed at lower d p can be explained by two factors: smaller particles result in lower terminal velocity, thus longer residence time within the reactor (see the left y-axis in Fig. 6(a)); second, smaller particles lead to faster reaction kinetics (see Eq. (28)), thus higher peak local

(a)

(c)

(b)

(d)

 3 = 1 104 mol m-2s-1 (left y Fig. 6. Effect of particle diameter on (a) average particle residence time at nsf,K axis) and average reduction extent at varying reaction rate constants (right y axis), and on centerline profiles of (b) oxygen partial pressure, (c) particle non-stoichiometry state, and (d) local reaction rate along z

 3 = 1 104 mol m-2s-1 for the baseline case; the thermodynamic prediction direction at the rate constant of nsf,K and local equilibrium conditions are added to facilitate our analysis.

reaction rates (see Fig. 6(d)). This also explains the trend observed in Fig. 6(a) that a higher reaction rate constant corresponds to a larger d p,cr .

4.3 Effect of selected operational parameters

(a)

(b)

(c) Fig. 7. Schematic of average reduction extent as a function of (a) particle mass flow rate at varying reaction rate constants, (b) gas mass flow rate at varying reaction rate constants, and (c) reaction rate constant at two cases of particle and gas flow rates while holding mg,in / mp,in unchanged; the thermodynamic prediction is included to facilitate our discussion.

Finally, the numerical model is employed to study the effects of selected operational parameters, particle and gas mass flow rates, on average outlet reduction extent. Investigating the effect of operating temperature does not yield much information and is therefore not discussed in the main text; we put this result in the Supplementary Information for interested readers. Fig. 7 shows the effects of particle and gas mass flow rates as well as reaction rate constant on average outlet reduction extent. For all cases, the maximum average particle volume fraction is found to be less

than 2  10 4 , confirming the assumption of dilute particle flow regime. Gas mass flow rate has a more pronounced effect on the reduction extent than particle mass flow rate at slow reaction kinetics; see Figs. 7(a) and (b). Another noteworthy finding is that when the mass flow rate ratio ( mg,in / mp,in ) is held constant, increase in absolute mass flow rates will lead to higher reduction extent in the chemical kinetics-limited region (see Fig. 7(c)). This is mainly due to the increased particle residence time arising from higher gas flow rates by effectively slowing down the particle velocity. In the gas advection-limited region, the reduction extent converges to the thermodynamic prediction as the reaction rate becomes higher, as shown in Figs. 7(a) and 7(b). Lower particle flow rates and/or higher gas flow rates will increase the gas advection capability, thus leading to higher thermodynamic reduction extents as seen in Figs. 7(a) and (b), respectively. Similar trends have been reported by Welte et al. [7] though under different operating conditions, suggesting relatively fast reaction kinetics in their experiments.  3 on  red To quantify the effects of mp,in , mg,in and nsf,K

p ave

as observed in Fig. 7, an empirical

correlation is developed as a function of three pertinent dimensionless numbers based on nonlinear regression using 197 numerical data points, each representing the outcome of a simulation:

 red

p ave,correlated

 1 2 0.1360   X 2.761  0.6523 X 2  12.68 X 2  219.0   0.05635 X 2  7.610  3   1.415   0.7484  X3  0.5634  0.1360  X  X 1 3

(36)

The three dimensionless numbers that appear are defined as:

X1 

mp,in mg,in

Reg 

X 2  Reg 

4mp,in

 Dt,in g

4mg,in

 Dt,in g

,

,

(37)

(38)

X 3  Da  Reg  3  nsf,K

M O2 Dt 2

g d p

1/2   hO02 ( ref )  Tp sO02 ( ref )   pO2 ,in   ,  exp        2 RT p p ref       3 ref

(39)

where Re g is the Reynolds number based on gas flow rate and Da is the Damköhler number calculated as the ratio of maximum reaction rate to gas phase advection. Note that the reference non-stoichiometry state appearing in Eq. (39) is specified as  ref  1105 in order to avoid the zero 3 value arising from the term  ref when  equals zero. The application range of this correlation is:

X1  (26.54,159.23) , X 2  (2.65,10.6) and X 3  (18.85, ) with other conditions held at their baseline values. The Reynolds number for the gas phase Re g (i.e., X2) is in the range of (2.65, 10.6), thus justifying the assumption (ii) of laminar gas flow as made in Section 2. A parity plot for the correlation versus the numerical results can be found from the Supplementary Information and reveals good agreement—R2=0.994. This correlation provides a means of quantifying the influences of gas and particle mass flow rates on the reduction extent for both slow and fast reaction kinetics circumstances without the need for detailed numerical simulations. 5

Conclusions A steady-state numerical model coupling mass and momentum transfer as well as reaction kinetics based on the Euler–Lagrange approach has been applied to an isothermal flow reactor composed of a downward flow of ceria particles undergoing reduction counter to an upward gas flow. The effects of reaction kinetics as well as design and operational parameters on reduction extent were investigated. The main factor that limits the reduction extent has been identified under varying conditions.

The reduction extent is found to increase with higher reaction rate constant until achieving the thermodynamic upper limit at a critical value. This critical rate constant is indicative of a turning point on the conversion-limiting factor from a chemical kinetics-limited regime to a regime that is limited by gas advection. Reactor design and operational parameters were also explored. The reactor length affects the reduction extent by influencing the particle residence time. There exists a critical reactor length above which the reduction extent always remains at the thermodynamic upper limit. The particles suffer from insufficient reaction time, and a lower reduction extent is achieved, when the reactor tube length is less than its critical value. Particle diameter influences reduction extent through both particle residence time and reaction kinetics. Smaller particles are characterized by longer residence time and faster chemical kinetics, so a critical particle diameter exists below which there is no further gain in the reduction extent. Caution must be taken to avoid particle entrainment for smaller particles at higher gas flow rates in order to realize CF configuration. An empirical correlation has been developed to quantify the influences of particle and gas mass flow rates at varying reaction kinetics. Future work will include the addition of a two-phase heat transfer model to predict more realistic temperature profiles as encountered in practical operation and to check/identify under what conditions the present isothermal assumption is a good treatment. Experimental work is also expected to proceed in parallel in order to validate the model results. Nomenclature

Aox

frequency factor for ceria oxidation (bar-n)

Ap

surface area for a single particle (m2)

Ared

frequency factor for ceria reduction (s-1)

At

inner cross-section area for a reactor tube (m2)

Cdrag

drag coefficient (–)

Ci

volume-based molar concentration of species i (mol m-3)

Da

Damköhler number (–)

DO2 ,eff

effective mass diffusion coefficient (m2 s-1)

DO2 ,g

oxygen mass diffusion coefficient in the gas phase (m2 s-1)

DO,p

ambipolar mass diffusion coefficient of oxygen within ceria (m2 s-1)

dp

particle diameter (m)

Dt,in

inner diameter of the tube reactor (m)

Ea,ox , Ea,red activation energy for ceria oxidation and reduction (kJ mol-1)

ew

coefficient of restitution for inelastic particle collision with the wall (–)

f drag

drag factor (–)

Fother,i

term representing other forces appearing in Eq. (8) (kg m s-2)

Fs

solid mass flux (kg m-2 s-1)

fv

volume fraction (–)

gi

acceleration of gravity in the i direction (m s-2)

  kox,Bulfin , kred,Bulfin reaction rate constant for ceria oxidation and reduction appearing in Eq. (26) (s-1

bar-n, s-1)  2 , kox,K  3 reaction rate constant for ceria oxidation appearing in Eqs. (24) and (25), respectively, kox,K

(mol-1 m4 s-1 atm-1/2, mol-2 m7 s-1 atm-1/2)  2 , kred,K  3 reaction rate constant for ceria reduction appearing in Eqs. (24) and (25), respectively, kred,K

(mol-1 m4 s-1, mol-2 m7 s-1)

Lt

reactor tube length (m)

m

mass flow rate (kg s-1)

M

molar weight (kg mol-1)

mp

mass for a single particle (kg)

n

a constant parameter appearing in Eqs. (26) and (29) (–)

ni

molar flow rate of material i (mol s-1)

nO 2

molar quantify of oxygen appearing in Eq. (4) (mol)

NBulfin

rate modification factor introduced in Eq. (29) (–)

N total

total number of particle stream (–)

Np

total number of particles in a computational cell (–)

N p,j

particle number flow rate of stream j (s-1)

 2 , nsf,K  3 rate constant appearing in Eq. (27) and Eq. (28), respectively (mol m-2 s-1) nsf,K p , pO 2

pressure, oxygen partial pressure (atm)

pO* 2

* 1 dimensionless oxygen partial pressure as defined by pO2  pO2 ( psys  pO2 ) (–)

r

arbitrary radial location (m)

R

universal gas constant (J mol-1 K-1)

Rij

volume-averaged Reynolds stress appearing in Eq. (21) (m2 s-2)

Rt,in

inner radius of the reactor tube (m)

rO2 , rO2

molar oxygen production rate for ceria reduction defined based on unit particle surface area and unit particle volume, respectively, (mol m2 s-1, mol m3 s-1)

Re g

Reynolds number calculated based on gas flow rate appearing in Eq. (38)

Re r

relative Reynolds number defined in Eq. (11)

S DPM,mass

mass source term appearing in Eqs. (18) and (19) due to interphase mass transfer (

kg m-3s-1 ) S DPM,momentum,i

momentum source term appearing in Eq. (21) due to interphase drag force as well

-2 -2 as mass transfer ( kg m s )

t

time (s)

T

temperature (K)

U

superficial velocity (m s-1)

Vc

volume of a computation cell (m3)

Vp

volume of a single particle (m3)

v j ,i

velocity of phase j in the i direction (m s-1)

v j ,i

velocity fluctuation of phase j in the i direction (m s-1)

x

a constant parameter appearing in Eq. (29) (–)

xi

spatial coordinate in the i direction (m)

X 1 , X 2 , X 3 dimensionless numbers appearing in Eq.(36) (–) Yi

mass fraction of species i (–)

Greek Symbols

 impact

particle impact angle with the wall (  )



particle non-stoichiometry state (–)

 ij

isotropic tensor Kronecker delta (–)



change in particle non-stoichiometry state (–)

hO2

standard enthalpy change for ceria oxidation (J mol-1)

nO 2

change of oxygen molar flow rate from gas inlet to outlet (mol s-1)

sO2

standard entropy change for ceria oxidation (J mol-1 K-1)

t j

time duration for particle stream j to traverse a computational cell (s)

t p

time step used to discretize time domain of particle phase (s)

p,cell

step length factor used for particle time discretization (–)



dynamic viscosity (kg m-1 s-1)

 g,bulk

coefficient of bulk viscosity for the gas phase (kg m-1 s-1)



density (kg m-3)

p

particle relaxation/response time (s)

Subscripts



gas mainstream properties

atm

atmospheric

ave

average

correlated predicted from the correlation

cr

critical

eff

effective

eq

equilibrium state

g

gas phase

in

inlet

norm

in the normal direction

out

outlet

ox

oxidation

p

particle phase

r

relative

red

reduction

ref

referenced state

s

properties at particle surface

sys

system

t

tube

tangent

in the tangential direction

ter

terminal

thermodyn thermodynamic prediction

Superscripts 0,1

before and after particle collision with a certain boundary

Other Symbols superficial average intrinsic average over phase

CeCe

electron localized on a Ce ion appearing in Eqs. (24) and (25)

CeCe

Ce ion on its normal site with neutral charge appearing in Eqs. (24) and (25)

OO

oxygen ion on its normal site with neutral charge appearing in Eqs. (24) and (25)

VO

general expression for an oxygen vacancy appearing in Eq. (26)

VO , VO

singly and doubly ionized oxygen vacancy appearing in Eqs. (24) and (25), respectively

Abbreviations

CF

countercurrent flow

DPM

discrete phase model

FVM

finite volume method

ODE

ordinary differential equation

PDE

partial differential equation

PF

parallel or cocurrent flow

SIMPLE Semi-Implicit Method for Pressure-Linked Equations Conflict of interest None declared. Acknowledgement The financial support of the China Scholarship Council (Sha Li, grant no. [2015]3022, 201506020092) and the Australian Research Council (Wojciech Lipiński, Future Fellowship, award no. FT140101213) is gratefully acknowledged. We thank Dr Peter Kreider for the fruitful discussions of chemical kinetics models. Supplementary Information Supplementary information to this article can be found online at website link. References [1]

Chueh, W. C., Falter, C., Abbott, M., Scipio, D., Furler, P., Haile, S. M., et al., 2010, "High-flux solar-driven thermochemical dissociation of CO2 and H2O using nonstoichiometric ceria." Science, 330(6012), pp. 1797–1801.

[2]

[3] [4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

Marxer, D., Furler, P., Takacs, M., Steinfeld, A., 2017, "Solar thermochemical splitting of CO2 into separate streams of CO and O2 with high selectivity, stability, conversion, and efficiency." Energy & Environmental Science, 10(5), pp. 1142–1149. Lapp, J., Davidson, J. H., Lipiński, W., 2012, "Efficiency of two-step solar thermochemical nonstoichiometric redox cycles with heat recovery." Energy, 37(1), pp. 591–600. Bader, R., Venstrom, L. J., Davidson, J. H., Lipiński, W., 2013, "Thermodynamic Analysis of Isothermal Redox Cycling of Ceria for Solar Fuel Production." Energy & Fuels, 27(9), pp. 5533– 5544. Tou, M., Michalsky, R., Steinfeld, A., 2017, "Solar-Driven Thermochemical Splitting of CO2 and In Situ Separation of CO and O2 across a Ceria Redox Membrane Reactor." Joule, 1(1), pp. 146– 154. Scheffe, J. R., Welte, M., Steinfeld, A., 2014, "Thermal Reduction of Ceria within an Aerosol Reactor for H2O and CO2 Splitting." Industrial & engineering chemistry research, 53(6), pp. 2175– 2182. Welte, M., Barhoumi, R., Zbinden, A., Scheffe, J. R., Steinfeld, A., 2016, "Experimental Demonstration of the Thermochemical Reduction of Ceria in a Solar Aerosol Reactor." Industrial & engineering chemistry research, 55(40), pp. 10618–10625. Lapp, J., Lipiński, W., 2014, "Transient Three-Dimensional Heat Transfer Model of a Solar Thermochemical Reactor for H2O and CO2 Splitting Via Nonstoichiometric Ceria Redox Cycling." Journal of Solar Energy Engineering, 136(3), pp. 031006. Li, S., Kreider, P. B., Wheeler, V. M., Lipiński, W., 2019, "Thermodynamic Analyses of Fuel Production Via Solar-Driven Ceria-Based Nonstoichiometric Redox Cycling: A Case Study of the Isothermal Membrane Reactor System." Journal of Solar Energy Engineering, 141(2), pp. 021012. Li, S., Wheeler, V. M., Kreider, P. B., Lipiński, W., 2018, "Thermodynamic Analyses of Fuel Production via Solar-Driven Non-stoichiometric Metal Oxide Redox Cycling. Part 1. Revisiting Flow and Equilibrium Assumptions." Energy & Fuels, 32(10), pp. 10838–10847. Li, S., Wheeler, V. M., Kreider, P. B., Bader, R., Lipiński, W., 2018, "Thermodynamic Analyses of Fuel Production via Solar-Driven Non-stoichiometric Metal Oxide Redox Cycling. Part 2. Impact of Solid–Gas Flow Configurations and Active Material Composition on System-Level Efficiency." Energy & Fuels, 32(10), pp. 10848–10863. Francis, T. M., Lichty, P. R., Weimer, A. W., 2010, "Manganese oxide dissociation kinetics for the Mn2O3 thermochemical water-splitting cycle. Part 1: Experimental." Chemical Engineering Science, 65(12), pp. 3709–3717. Loutzenhiser, P. G., Elena Gálvez, M., Hischier, I., Graf, A., Steinfeld, A., 2010, "CO2 splitting in an aerosol flow reactor via the two-step Zn/ZnO solar thermochemical cycle." Chemical Engineering Science, 65(5), pp. 1855–1864. Z’Graggen, A., Steinfeld, A., 2009, "Heat and mass transfer analysis of a suspension of reacting particles subjected to concentrated solar radiation—Application to the steam-gasification of carbonaceous materials." International Journal of Heat and Mass Transfer, 52(1–2), pp. 385–395. Bader, R., Gampp, L., Breuillé, T., Haussener, S., Steinfeld, A., Lipiński, W., 2018, "Unsteady Radiative Heat Transfer Model of a Ceria Particle Suspension Undergoing Solar Thermochemical Reduction." Journal of Thermophysics and Heat Transfer, 33(1), pp. 63–77. Martinek, J., Bingham, C., Weimer, A. W., 2012, "Computational modeling of a multiple tube solar reactor with specularly reflective cavity walls. Part 2: Steam gasification of carbon." Chemical Engineering Science, 81, pp. 285–297. Haussener, S., Hirsch, D., Perkins, C., Weimer, A., Lewandowski, A., Steinfeld, A., 2009, "Modeling of a Multitube High-Temperature Solar Thermochemical Reactor for Hydrogen Production." Journal of Solar Energy Engineering, 131(2), pp. 024503. Groehn, A. J., Lewandowski, A., Yang, R., Weimer, A. W., 2016, "Hybrid radiation modeling for multi-phase solar-thermal reactor systems operated at high-temperature." Solar Energy, 140, pp. 130–140.

[19] [20]

[21] [22]

[23]

[24]

[25] [26]

[27] [28] [29] [30] [31] [32] [33] [34] [35]

[36] [37] [38] [39] [40]

[41] [42]

Lipiński, W., Steinfeld, A., 2005, "Transient radiative heat transfer within a suspension of coal particles undergoing steam gasification." Heat and Mass Transfer, 41(11), pp. 1021–1032. Lipiński, W., Thommen, D., Steinfeld, A., 2006, "Unsteady radiative heat transfer within a suspension of ZnO particles undergoing thermal dissociation." Chemical Engineering Science, 61(21), pp. 7029–7035. Maag, G., Lipiński, W., Steinfeld, A., 2009, "Particle–gas reacting flow under concentrated solar irradiation." International Journal of Heat and Mass Transfer, 52(21–22), pp. 4997–5004. Lipiński, W., Davidson, J. H., Haussener, S., Klausner, J. F., Mehdizadeh, A. M., Petrasch, J., et al., 2013, "Review of Heat Transfer Research for Solar Thermochemical Applications." Journal of Thermal Science and Engineering Applications, 5(2), pp. 021005. Perkins, C., Weimer, A., 2007, "Computational Fluid Dynamics Simulation of a Tubular Aerosol Reactor for Solar Thermal ZnO Decomposition." Journal of Solar Energy Engineering, 129(4), pp. 391–404. Francis, T. M., Perkins, C., Weimer, A. W., 2010, "Manganese oxide dissociation kinetics for the Mn2O3 thermochemical water-splitting cycle. Part 2: CFD model." Chemical Engineering Science, 65(15), pp. 4397–4410. Wheeler, V. M., Bader, R., Kreider, P. B., Hangi, M., Haussener, S., Lipiński, W., 2017, "Modelling of solar thermochemical reaction systems." Solar Energy, 156, pp. 149–168. Abanades, S., Charvin, P., Flamant, G., 2007, "Design and simulation of a solar chemical reactor for the thermal reduction of metal oxides: Case study of zinc oxide dissociation." Chemical Engineering Science, 62(22), pp. 6323–6333. Oles, A. S., Jackson, G. S., 2015, "Modeling of a concentrated-solar, falling-particle receiver for ceria reduction." Solar Energy, 122, pp. 126–147. Crowe, C., Schwarzkopf, J., Sommerfeld, M., Tsuji, Y., Multiphase Flows with Droplets and Particles Boca Raton: CRC Press, 2012. Schiller, L., Nauman, A., 1933, "On the Basic Computations in Gravity Processing." Z Vereines Deutscher Inge, 77, pp. 318–321. Clift, R., 1970, "The Motion of Particles in Turbulent Gas-Streams." Proc Chemeca'70, 1, pp. 14. Putnam, A., 1961, "Integratable form of droplet drag coefficient." 31(10), pp. 1467–1468. Morsi, S., Alexander, A., 1972, "An investigation of particle trajectories in two-phase flow systems." Journal of Fluid mechanics, 55(2), pp. 193–208. Whitaker, S., The method of volume averaging. Springer Science & Business Media, 2013. Bird, R. B., Stewart, W. E., Lightfoot, E. N., Klingenberg, D. J., Introductory transport phenomena. Wiley Global Education, 2015. Ackermann, S., Scheffe, J. R., Steinfeld, A., 2014, "Diffusion of Oxygen in Ceria at Elevated Temperatures and Its Application to H2O/CO2 Splitting Thermochemical Redox Cycles." The Journal of Physical Chemistry C, 118(10), pp. 5216–5225. Kundu, P. K., Dowling, D. R., Tryggvason, G., Cohen, I. M., Fluid mechanics. 2015. Graham, T., 1846, "On the motion of gases." Philosophical Transactions of the Royal Society of London, (136), pp. 573–631. Davidson, T. A., 1993, "Simple and accurate method for calculating viscosity of gaseous mixtures(Report of Investigations, 1993)." pp. 1–12. Yaws, C. L., Transport properties of chemicals and hydrocarbons. William Andrew, 2014. Keene, D. J., Davidson, J. H., Lipiński, W., 2013, "A Model of Transient Heat and Mass Transfer in a Heterogeneous Medium of Ceria Undergoing Nonstoichiometric Reduction." Journal of Heat Transfer, 135(5), pp. 052701. Keene, D. J., Lipiński, W., Davidson, J. H., 2014, "The effects of morphology on the thermal reduction of nonstoichiometric ceria." Chemical Engineering Science, 111, pp. 231–243. Bulfin, B., Lowe, A. J., Keogh, K. A., Murphy, B. E., Lübben, O., Krasnikov, S. A., et al., 2013, "Analytical Model of CeO2 Oxidation and Reduction." The Journal of Physical Chemistry C, 117(46), pp. 24129–24137.

[43] [44]

[45]

[46] [47] [48]

[49] [50] [51] [52]

Panlener, R., Blumenthal, R., Garnier, J., 1975, "A thermodynamic study of nonstoichiometric cerium dioxide." Journal of Physics and Chemistry of Solids, 36(11), pp. 1213–1222. Bulfin, B., Hoffmann, L., de Oliveira, L., Knoblauch, N., Call, F., Roeb, M., et al., 2016, "Statistical thermodynamics of non-stoichiometric ceria and ceria zirconia solid solutions." Physical chemistry chemical physics : PCCP, 18(33), pp. 23147–23154. Chueh, W. C., Haile, S. M., 2010, "A thermochemical study of ceria: exploiting an old material for new modes of energy conversion and CO2 mitigation." Philosophical transactions Series A, Mathematical, physical, and engineering sciences, 368(1923), pp. 3269–3294. Laín, S., Sommerfeld, M., 2008, "Euler/Lagrange computations of pneumatic conveying in a horizontal channel with different wall roughness." Powder Technology, 184(1), pp. 76–88. ANSYS, Inc., ANSYS Fluent Theory Guide 17.1. Hernandez-Perez, V., Abdulkadir, M., Azzopardi, B., 2011, "Grid generation issues in the CFD modelling of two-phase flow in a pipe." The Journal of Computational Multiphase Flows, 3(1), pp. 13–26. Laín, S., Sommerfeld, M., 2013, "Characterisation of pneumatic conveying systems using the Euler/Lagrange approach." Powder Technology, 235, pp. 764–782. Patankar, S., Numerical heat transfer and fluid flow. CRC press, 2018. Belloni, A., Industrial gases processing. John Wiley & Sons, 2008. Brendelberger, S., Roeb, M., Lange, M., Sattler, C., 2015, "Counter flow sweep gas demand for the ceria redox cycle." Solar Energy, 122, pp. 1011–1022.

Numerical modelling of ceria undergoing reduction in a particle–gas counter-flow: effects of chemical kinetics under isothermal conditions Sha Lia, Vincent M. Wheelerb, Apurv Kumara, c, and Wojciech Lipińskia,* a

Research School of Electrical, Energy and Materials Engineering, The Australian National University, Canberra,

ACT 2601, Australia b

Department of Engineering and Technology, University of Wisconsin–Stout, Menomonie, WI 54751, USA

c

School of Science, Engineering and Information Technology, Federation University, Mt. Helen, Ballarat, Australia

Highlights 

Particle–gas reactive counter-flow analyzed using Euler–Lagrange approach;



Factors limiting the outlet reduction extent identified at varying conditions;



Design strategies proposed to approach the upper limit of reaction extent;



Correlation developed relating reduction extent to flow rate and kinetics.

Numerical modelling of ceria undergoing reduction in a particle–gas counter-flow: effects of chemical kinetics under isothermal conditions Sha Lia, Vincent M. Wheelerb, Apurv Kumara, c, and Wojciech Lipińskia,* a

Research School of Electrical, Energy and Materials Engineering, The Australian National University, Canberra,

ACT 2601, Australia b

Department of Engineering and Technology, University of Wisconsin–Stout, Menomonie, WI 54751, USA

c

School of Science, Engineering and Information Technology, Federation University, Mt. Helen, Ballarat, Australia

Credit author statement Sha Li:

Conceptualization, Methodology, Investigation, Writing – Original Draft, Funding acquisition

Vincent M. Wheeler: Conceptualization, Writing - Review & Editing, Supervision Apurv Kumar:

Conceptualization, Writing – Review & Editing, Supervision

Wojciech Lipiński:

Conceptualization, Writing – Review & Editing, Supervision, Project administration, Funding acquisition

We declare no conflict of interest.

1