Modeling of Ceria Reduction in a Solar Thermochemical Reactor via DEM Method
Han Zhang , Joseph D. Smith ^{*}
Department of Chemical and Biochemical Engineering, Missouri University of Science and Technology, Rolla, MO 65409, United States
* Correspondence: Joseph D. Smith
Academic Editor: Joaquin AlonsoMontesinos
Special Issue: Photovoltaic Solar Systems and Solar Thermal Plants
Received: November 09, 2020  Accepted: January 19, 2021  Published: January 30, 2021
Journal of Energy and Power Technology 2021, Volume 3, Issue 1, doi:10.21926/jept.2101008
Recommended citation: Zhang H, Smith JD. Modeling of Ceria Reduction in a Solar Thermochemical Reactor via DEM Method. Journal of Energy and Power Technology 2021; 3(1): 008; doi:10.21926/jept.2101008.
© 2021 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.
Abstract
Solar thermochemical reactor provides an attractive approach that utilizes the most common solar radiation as the thermal driving force to motivate the reaction between CO_{2} and metal oxides, which is also called metal oxide redox pairbased thermochemical cycles. The CeO_{2}/CeO_{2δ} is widely used in the twostep redox process due to its advantages including fastredox kinetics, high crystallographic stability of a wide range of reacting oxygen nonstoichiometry, and relatively high oxygen solidstate conductivity. In this work, a threedimensional transient numerical analysis has been completed to study the performance of a CeO_{2} reduction reaction in a 1/8^{th} segment region of a novel partition cavityreceiver reactor. The porous CeO_{2} catalyst was analyzed using the discrete element method (DEM) to capture the heat transfer and reactive performances. The catalyst textural properties (particle size and void fraction) and reaction conditions (gas flow rate and radiative power input) were investigated in the CeO_{2} reduction reaction. The results indicated that increasing the catalyst specific surface area and the temperature are beneficial to the O_{2} production and further CO_{2} conversion.
Graphical abstract
Keywords
Solar thermochemical reactor; ceria reduction; radiation; transient; DEM; chemical energy storage
1. Introduction
With the growing energy demands in the world, environmental protection issues have attracted a lot of attention related to energy technology development. Among various energy resources, fossil fuel remains the primary and stable energy source. Along with the use of fossil fuels, carbon dioxide emission rate is a serious factor that affects climate change and the global environment. Figure 1 shows global carbon dioxide emissions from 1990 to 2019 [1].
Figure 1 Global carbon dioxide emissions from 1990 to 2019 [1].
From this figure, CO_{2} emission continues at a high level in the recent three years. However, the emission of 2019 kept flat in 2019 due to the expanding application of renewable energy. Thanks to the numerous studies focusing on using clean energies or novel technologies, they showed important effect on the reduction of CO_{2} emission [2,3,4,5,6]. Since CO_{2} is also a carbon source and can be converted into useful and variable materials, work has been focused on finding efficient technology to accomplish this process. Decades ago, people began studying water electrolysis to produce hydrogen, since water is a cheap and available resource. However, hydrogen storage is an issue due to its high volatile property [7]. Compared with water electrolysis, syngas production from CO_{2} and H_{2}O does not only provide an option for environmental protection and transportation safety but also offers an attractive energy fuel for further applications during the chemical process, like FischerTropsch technology [8]. Solar energy is essentially a green technology with inexhaustible power production potential. Under the concentrated energy from the sun, CO_{2} and H_{2}O can be converted to syngas using proper metal oxide redox pairs such as ZnO/Zn [9], Fe_{3}O_{4}/FeO [10], and many more advanced couples [11]. The appearance of producing syngas via solar energy provides a brandnew idea for both solar energy application and chemical fuel production.
The twostep solar thermochemical process is a redox cycle based on metal oxide redox. This chemical process is represented as the hightemperature reduction (see Eq. (1)) and lowtemperature oxidation (see Eq. (2) and (3)).
Hightemperature reduction:
\[ M_xO_y\ \rightarrow\ \ M_xO_{y\delta}+\ \frac{\delta}{2}\ O_2\tag{1} \]
Lowtemperature oxidation:
\[ M_xO_{y\delta}+\ \delta{CO}_2\ \rightarrow\ M_xO_y+\ \delta CO\tag{2} \]
\[ M_xO_{y\delta}+\ \delta H_2O\ \rightarrow\ M_xO_y+\ \delta H_2\tag{3} \]
In the first step, the metal oxide is reduced to form pure metal or reductive metal oxide that releases oxygen over 1273 K [12]. This reaction is endothermic, so excess energy must be supplied using concentrated solar energy. In the second step, the oxidation reactions are exothermic and occur at a lower temperature (<1000 K). The feed CO_{2}/H_{2}O react with metal or reduced metal oxide and lose oxygen to form CO/H_{2}. Correspondingly, the metal or reduced metal oxide is oxidized to the former oxidation state. Compared with syngas production cycles, the most remarkable advantage of the twostep solar thermochemical process is that the water splitting and carbon dioxide splitting share the same reaction mechanism, which provides a great convenience for further analysis [13]. The reaction formula are shown in Eq. (2) and Eq. (3).
The metal oxide pairs for the twostep solar thermochemical redox process can be divided into two groups, volatile and nonvolatile. The classification is based on the material physical properties and whether phase change is occurred during the metal oxide reduction [14]. Since the temperature requirement for reduction is in the range of 1500 K to 2000 K and up to 2300 K, possibly higher than the boiling temperature of the reduced species, the volatile metal oxide may experience a phase change. There are some representative metal oxide pairs, including ZnO/Zn, CdO/Cd, SnO_{2}/SnO, GeO_{2}/GeO, and so forth [15]. Even though the volatile redox pair results in an extra entropy gain for the system due to the phase transition, the recombination of oxygen and reduced metal or metal oxide decreases the available reduced metal or metal oxide in the oxidization process, which has a negative impact on fuel production of volatile pairs [16]. Nonvolatile redox pairs (e.g. ferrites, ceria, hercynite, and perovskites) [17,18] can maintain their solid state at high temperatures, which avoid the recombination issue. Other morphology characteristics of nonvolatile redox pairs, such as pore size, porosity, and specific surface area, also affect the reaction rate.
Cerium oxide is a representative redox catalyst used to produce carbon monoxide and hydrogen in the early 1980s due to its excellent performance in releasing and recovering oxygen [19]. The application of the nonstoichiometric CeO_{2} redox pair for solar thermochemical reactions flourished in the last decade. In 2006, Flamant et al. [20] proposed the application of CeO_{2}/Ce_{2}O_{3} redox pair to split water. The reduction of CeO_{2} was proceeded at 2000 °C, 100~200 mbar, followed by the oxidation step in the temperature range of 400600 °C. Since reduction temperature is high enough to transfer CeO_{2} to a molten state, which strongly reduces the catalyst activity and quality, the study of CeO_{2} was focused on how to avoid the phase transfer under a high reduction temperature. Nonstoichiometric CeO_{2} (CeO_{2δ}) was proposed to solve the issue of melting at high temperature. Based on the twostep thermochemical process, the reactions can be represented by
\[ CeO_2 \to CeO_{2\delta}+ \frac{\delta}{2} O_2\tag{4} \]
\[ CeO_{2\delta}+ \delta CO_2 \to CeO_2+ \delta CO\tag{5} \]
\[ CeO_{2\delta}+ \delta H_2 O \to CeO_2+ \delta H_2\tag{6} \]
The net reaction enthalpies of CO_{2} splitting ($ CO_2\to CO+\frac{1}{2}O_2 $) and H_{2}O splitting ($ H_2 O\to H_2+\frac{1}{2}O_2 $) are 283 kJ per mole of CO_{2} and 242 kJ per mole of H_{2}O, respectively [21,22]. At 1500 and 0.1 mbar O_{2} partial pressure, the enthalpy change of ceria reduction is around 475 kJ per half mole of O_{2} [21]. Therefore, the reaction enthalpies to generate one mole of CO and H_{2} are estimated 192 kJ/mol and 233 kJ/mol.
The solartofuel efficiency (see Eq. (7)) is used to show the reaction performance [22].
\[ \eta_{solartofuel}=\frac{\eta_{fuel}HHV}{Q_{solar}}\tag{7} \]
Chueh et al. [13] presented the application of porous monolithic ceria on the twostep solar thermochemical CO_{2}/H_{2}O splitting. The results showed that porous monolith ceria can reach a stable fuel generation state after 500 redox cycles with an estimated 0.8% peak solartofuel efficiency. Furler et al. [23] proposed using porous ceria felt as the catalyst to support solar CO_{2}/H_{2}O splitting. The O_{2} and fuel production are 2.89±0.27 mL g1 CeO_{2} and 5.88±0.43 mL g1 CeO_{2}. Furler et al. [24] also proposed a reticulated porous ceria (RPC) foam structure for redox reactions, which derived a 3.53% peak solartofuel efficiency. The efficiency improved more than four times compared with former porous monolithic ceria. Venstrom et al. [25] presented a novel morphology of CeO_{2}, threedimensionalordered macroporous (3DOM) CeO_{2}, which has both higher specific surface area and effective pore interconnections. The fuel production rate was enhanced by as much as 260% and 175% compared with traditional low porosity ceria and nonordered mesoporous ceria. Oliveira et al. [26] further studied the 3DOM ceria structure for solar thermochemical CO_{2} splitting in a tubular reactor. The result showed that the maximum fuel production rate is three times higher than the dualscale porosity ceria reticulated porous foam.
Based on the difference of heat transfer type, solar thermochemical reactors can be classified as indirectlyirradiated reactors and directlyirradiated reactors. For indirectlyirradiated reactors, the reactor chamber material is generally made from opaque metals. Incoming solar radiant flux to the external wall of the reactor is transferred to heat conduction to raise the temperature of the metal oxide. For directlyirradiated reactors, the solar radiation can directly touch the surface of the metal oxide and use the incident radiation to heat up the oxide. Since the radiation term in the energy conservation equation also participates in heat transfer in directlyirradiated reactors, the reactor can achieve a higher temperature. To fully utilize the solar radiation, the cavityreceiver reactor is commonly used for directlyirradiated solar thermochemical processes. In a cavityreceiver reactor, the space provided in the cavity provides or enhances radiation to the catalyst inner surface, which forms an approximate blackbody. Different types of cavityreceiver reactors have been discussed and used over the past several years to prove that direct irradiation has a huge advantage on the twostep solar H_{2}O/CO_{2} splitting process [27].
To better predict and validate the CeO_{2} supported solar thermochemical twostep redox cycle, several chemical kinetics were studied and utilized to model the solar thermochemical process performances in recent years. Thermodynamic equilibrium model is a representative reaction model that widely used in the simulations of CeO_{2} twostep redox reactions. Calle et al. [28], Groehn et al. [29], Bulfin et al. [30] introduced the relation between nonstoichiometric coefficient with oxygen partial pressure and temperature [31] in their process simulation and reactor modeling, respectively. It was proved that the model has a higher agreement with the model of Zinkevich et al. [32]. Zoller et al. [26] used the thermodynamic equilibrium model to simulate a 50kW solar receiverreactor, which obtained the solartofuel efficiency larger than 10% with dualscale reticulated porous ceramic. The most basic chemical kinetics reaction model is based on the Arrhenius equation proposed by Ishida et al. [33]. This model was applied in different tubular reactor studies [34,35]. However, the accuracy needs to be improved, due to the limitation of the current reported values of coefficients in the Arrhenius equation [36]. Arifin et al. [37] and McDaniel et al. [38] used the improved solidstate kinetic model, which includes the influences of reaction order of the oxidant mole fraction and the progress of the rate of oxidation [39], With their studies, the solidstate kinetic model is further parametrized for accurately evaluating the reaction rates. Based on the theory of Kröger et al. [40], the crystal site relations of CeO_{2} reduction was discussed in the work of Scheffe et al. [41] and corresponding reaction mechanism was applied in the work of Bader et al. [42] to study the CeO_{2} reduction performances in a particle suspension system. Haeussler et al. [43,44] studied the twostep redox cycling of reticulated porous ceria structures in a monolithic cavitytype reactor. By increasing the reduction temperature or decreasing the operating pressure, the fuel production yields can reach up to 341 µmol/g and the peak solartofuel efficiency achieves up to 8%. Parthasarathy et al. [45] used CFD to validate a solar thermochemical reactor with ceria RPC structure. In the studied power range (0.8 kW to 3.8 kW), the temperature of the CFD model has an error range between 0.3% to 7.15% compared with experiment results. Kyrimis et al. [46] numerical studied the scaleup solar thermochemical reactor geometry. The results showed that the scaleup geometry has a benefit for the solartofuel efficiency at the beginning of the reduction but limited due to the thickness of CeO_{2} RPC.
Gassolid performance plays a vital role in the numerical simulation of solar thermochemical processes. The EulerianLagrangian (CFDDEM) method, which uses the DEM method to solve the transfer performances of particles individually in the gas continuum phase, is a representative approach to capture the variations of the dispersed phase. It was widely used in the simulations of fluidized bed reactors [47,48], however, the application in the solar thermochemical processes for hightemperature chemical reactions is still in the infantile stage. Bellan et al. did a series of work on the heat transfer performances in different types of fluidized bed reactors for solar gasification and solar thermochemical cycles using the CFDDEM method and the results showed high agreement between the experimental and simulation results [49,50]. Morris et al. used CFDDEM method to study the heat transfer performances in the novel solar receiver designed by the National Renewable Energy Laboratory and proved that the improved continuum model has a high accuracy to predict the heat transfer coefficient [51]. Although some studies have been reported to focus on the heat transfer performances in the solar fluidized bed reactors using the CFDDEM method, very few works consider the reactions in the solar thermochemical twostep metal redox cycles. Additionally, the transfer performances in the porous catalyst simulated by the DEM methods are worthy of discussing intensively.
In this paper, a threedimensional (3D) transient multiphysics computational fluid dynamics model is applied to study heat and mass transfer of a hightemperature CeO_{2} reduction process in a novel partition cavity solar thermochemical reactor. The model uses the discrete element method (DEM) [52] to approximate the structure of a CeO_{2} layer as a packed bed. Several parameters defining the operating conditions and the catalyst layer formation have been studied to investigate the performance of the CeO_{2} reduction in the novel partitioncavity thermochemical reactor.
2. Reactor Design Concept
The reactor (shown in Figure 2) consists of a cylindrical cavity and a solar radiation receiver. An annular cylindrical catalyst region (60 mm i.d., 100 mm o.d., and 150 mm height) is integrated into the cavity section. The reactive region is separated from the cavity center by a transparent quartz partition (145 mm height). A 5 mm gap between the bottom of the transparent partition and cavity bottom leaves space for fluid to pass through the entire length of the reactive packed bed. Additionally, the transparent quartz allows solar radiation to enter the reactive region. With the partition design, the reactive cavity not only provides a relatively long path for contact of gas and solid but also enhance radiation heat transfer to the reactive region. The reactive region is surrounded by a ceramic insulation material inside the stainlesssteel outer shell. Concentrated solar radiation enters the reactor through the 240 mm receiver window and focuses on the interface of the solar receiver and reactive cavity. Concentrated radiation distributes over the inner cavity surface and enters the reactive region. Eight inlets and outlets, each with a diameter of 8 mm, are evenly distributed around the top edge of the solar receiver and reactive cavity side walls (see Figure 2 and Figure 3).
Figure 2 Partition cavity solar reactor schematic diagram.
Figure 3 Simulation reactor model. (a) top view; (b) crosssection.
3. Approach
3.1 Reactor Model
Figure 3 shows the solar reactor modeling geometry in top view and crosssection. Since the reactor has a symmetric structure, only oneeighth of the entire geometry is simulated, with the assumption that radiation distribution and thermal losses are symmetric. Plane MN and MP are periodic boundary conditions with a 45o rotational angle. The stainlesssteel housing and receiver are considered insignificant resistances, which are included in the wall boundary conditions to approximate heat loss from the reactor. Cerium dioxide (CeO_{2}) is packed in the reactive region to form a porous media structure approximated by the DEM method. By controlling the DEM injecting conditions, the catalyst particles are uniformly distributed through the entire reactive region that creates the appropriate gas flow paths to effectively mimic the structure of a porous media region. The reactive section is housed in the ceramic insulation material Al_{2}O_{3}SiO_{2} with surfaces to exchange radiation heat transfer. Radiative heat transfer is implemented using the discrete ordinates model (DOM) [52] with gray thermal radiation model, which considers the radiative properties of participated media and the surrounding surfaces are identical for all wavelengths [53]. CeO_{2} is assumed to be a black body radiator and Al_{2}O_{3}SiO_{2} layer is assumed to be a good thermal insulator with an emittance of 0.28 [54]. The receiver side wall is highly reflective material with a reflectance of 0.9 [55]. The incident solar power is assumed to be evenly distributed on the partition surface for simplification. Argon (Ar), entering the reactor from radial inlets on the solar receiver, works as an inert purge gas during preheating and reduction processes to improve the heat convection and accelerates the motion of generated oxygen.
3.2 Numerical Analysis and Validation
The simulated reactor schematic diagram in Figure 3 (b) represents the computational domains considered in the reactor model. Generally, the domains can be divided into nonreactive regions and reactive regions. Nonreactive domains include the solar receiver, the reactor cavity, and the inlets and outlets, which are simulated as fluid regions with laminar flow and ideal gas mixture. The transparent quartz partition and Al_{2}O_{3}SiO_{2} insulation are considered as solid walls with corresponding thermal properties (conductivity and emittance). The reactive domain is a multiphase region gas mixture flowing through a packed bed of CeO_{2} particles. To set up the link between the particles and fluid phases, a twoway coupling method is introduced that solves the heat and mass transfer conservation equations in both phases [52].
3.2.1 Nonreactive Regions
The continuity, momentum, energy, and mass conservation equations of fluid regions areas follow:
\[ \frac{\partial\rho_f}{\partial t}+\nabla\cdot(\rho_f v)=0\tag{8} \]
\[ \frac{\partial(\rho_f v)}{\partial t}+\nabla\cdot(\rho_fvv)=\nabla p+\nabla\cdot\mu\left(\nabla{v}+(\nabla{v})^T\frac{2}{3}{I}\nabla\cdot{v}\right)+{f}_{b}\tag{9} \]
\[ \frac{\partial\left(\rho_fH_f\right)}{\partial t}+\nabla\cdot\left(\rho_fH_f{v}\right)=\nabla\cdot\left(k\nabla T\right)\tag{10} \]
\[ \frac{\partial\left(\rho_fY_i\right)}{\partial t}+\mathrm{\nabla}\cdot\left(\rho_f{v}Y_i\right)=0\tag{11} \]
where ρ_{f} is fluid density, v is velocity vector, p is pressure, μ is fluid dynamic viscosity,I is identity matrix, f_{b} is momentum body force (gravity), H_{f} is fluid enthalpy, k is fluid thermal conductivity, and Y_{i} is component concentration [52].
Since there is no flow in the solid regions, only the energy conservation equation is considered. Eq. (12) and (13) show the energy conservation of partition and insulation:
\[ \frac{\partial\left(\rho_{part}H_{part}\right)}{\partial t}+\mathrm{\nabla}\cdot\left(\rho_{part}H_{part}{v}\right)=\mathrm{\nabla}\cdot\left(k_{part}\mathrm{\nabla T}\right)+S_{E,rad}\tag{12} \]
\[ \frac{\partial\left(\rho_sH_s\right)}{\partial t}+\mathrm{\nabla}\cdot\left(\rho_sH_sv\right)=\mathrm{\nabla}\cdot\left(k_s\mathrm{\nabla T}\right)+\mathrm{\nabla}q_r\tag{13} \]
where ρ_{part} is partition density, ρ_{s} is insulation density, H_{part} is partition enthalpy, H_{s} is insulation enthalpy, k_{part} is partition thermal conductivity, k_{s} is insulation thermal conductivity, S_{E,rad} is source term of radiation, ∇q_{r} is reradiation term. The input radiation and the reradiation are included in S_{E,rad}, which can be solved from the radiative transfer equation:
\[ \frac{dI}{ds}=\beta I+\kappa_aI_b+\frac{\kappa_s}{4\pi}\int_{4\pi} I\left(\Omega\right)d\left(\Omega\right)+\kappa_{pa}I_{pb}+\frac{\kappa_{ps}}{4\pi}\int_{4\pi} I\left(\Omega\right)d\left(\Omega\right)\tag{14} \]
where I is radiative intensity, κ_{a} is absorption coefficient, κ_{s} is scattering coefficient, $I_b$ is black body intensity, Ω is solid angle, κ_{pa} is particle absorption coefficient, κ_{ps} is particle scattering coefficient,$I_{pb}$ is particle black body intensity, β is extinction coefficient, and s is distance in the Ω direction. As the Ar and O_{2} are low radiative materials, the radiative properties κ_{a} and κ_{s} are assumed to equal zero [56]. Reradiation is used to couple the fluid and radiant energy field, which is shown as:
\[ \mathrm{\nabla}q_r=\kappa_{pa}\left(4\pi I_b\int_{4\pi} I d\Omega\right)\tag{15} \]
3.2.2 Reactive Region
The multiphase domain is the corereactive region in the reactor consisting of two phases: gas mixture fluid phase and CeO_{2} particle Lagrangian phase. The fluid continuous phase should consider the interaction between gas mixture and CeO_{2} particles, which can be conveyed as source terms in governing equations. The conservation equations are as follows:
\[ \frac{\partial{(\alpha\rho}_f)}{\partial t}+\nabla\cdot(\alpha\rho_f{v})=S_{m,O_2}\tag{16} \]
\[ \frac{\partial({\alpha\rho}_f{v})}{\partial t}+\nabla\cdot({\alpha\rho}_f{vv})=\nabla\alpha p+\nabla\cdot\alpha\mu\left(\nabla{v}+(\nabla{v})^T\frac{2}{3}{I}\nabla\cdot{v}\right)+{\alpha{f}}_{b}+S_V\tag{17} \]
\[ \frac{\partial\left(\alpha\rho_fH_f\right)}{\partial t}+\mathrm{\nabla}\cdot\left({\alpha\rho}_fH_f{v}\right)=\mathrm{\nabla}\cdot\left(\alpha k\nabla T\right)+h_{fs}A_{fs}\left(T_sT_f\right)+\mathrm{\nabla}q_r{+S}_E\tag{18} \]
\[ \frac{\partial\left(\alpha\rho_fY_{O_2}\right)}{\partial t}+\mathrm{\nabla}\cdot\left({\alpha\rho}_f{v}Y_{O_2}\right)=S_{m,\ O_2}\tag{19} \]
where α is the void fraction of fluid phase,$S_{m, O_2}$, S_{V}, and S_{E} are the binary phase coupling mass source term, momentum source term, and energy source term, respectively. These terms can be expanded and shown as:
\[ S_V=\frac{1}{\Delta t}\sum_{\pi}{(\int_{t}^{t+\Delta t}\int_{V_c}{\delta\left({r}{r}_{\pi}\right)n_\pi\left({F}_{s}+{\dot{m}}_p{v}_{p}\right)dVdt)}}\tag{20} \]
\[ S_E=\frac{1}{\Delta t}\sum_{\pi}{(\int_{t}^{t+\Delta t}\int_{V_c}{\delta\left({r}{r}_{\pi}\right)n_\pi\left(Q_t+{F}_{s}\cdot v_p+\frac{1}{2}{\dot{m}}_p{{v}_{p}}^2+{\dot{m}}_ph\right)dVdt)}}\tag{21} \]
\[ S_{m,O_2}=\frac{1}{\Delta t}\sum_{\pi}{(\int_{t}^{t+\Delta t}\int_{V_c}{\delta\left({r}{r}_{\pi}\right)n_\pi{\dot{m}}_pdVdt)}}\tag{22} \]
where F_{s} is particle surface force and m_{p} is particle mass transfer rate. The volume integral is based on cells. The Dirac delta function accounts for the parcels not represented in the cell. The concrete formations of F_{s} and m_{p} are discussed with the governing equations for DEM particles. By taking the particle motion, energy, and mass influences into consideration, the fluid governing equations approximate fluid/particle interactions using the DEM packing conditions.
Since CeO_{2} particles work as discrete elements, the momentum conservation is shown based on surface and body force. Particle mass and heat transfer happen on interphase, which should be considered together. The governing equations are represented below:
\[ m_p\frac{d{v}_{p}}{dt}={F}_{s}+{F}_{b}\tag{23} \]
\[ m_pc_p\frac{dT_p}{dt}=Q_{conv}+Q_{rad}+Q_{reac}\tag{24} \]
\[ \frac{dm_p}{dt}={\dot{m}}_p\tag{25} \]
where F_{s} is surface force including drag force and pressure gradient, F_{b} is body force including gravity and contact force, Q_{conv} is convective heat transfer from the fluid phase to the particle, Q_{rad} is radiative heat transfer rate, Q_{reac} is heat of reduction reaction. Eq. (26) shows how Q_{conv} can be formulated:
\[ Q_{conv}=h_{fp}A_s\left(TT_p\right)\tag{26} \]
where h_{fp} is heat transfer coefficient, A_{s} is particle surface area, T is fluid temperature, T_{p} is particle temperature. Heat transfer coefficient h_{fp} can be derived from the relation of particle Nusselt number [57]:
\[ {Nu}_p=\frac{h_{fp}D_p}{k}\tag{27} \]
where k is particle thermal conductivity and D_{p} is particle diameter. Since Re_{p}<5000, the RanzMarshall correlation [58] can be applied to determine the particle Nusselt number:
\[ {Nu}_p=2\left(1+0.3{Re}_p^{{1}/{2}}{Pr}^{{1}/{3}}\right)\tag{28} \]
where Pr is from the fluid phase.
The particle participated radiative heat transfer Q_{rad} is defined as:
\[ Q_{rad}=\frac{A_s}{4}Q_{a,p}\left(G4\sigma T_p^4\right)\tag{29} \]
where Q_{a,p} is absorption efficiency of particle, G is incident radiative heat flux, and σ is the StefanBoltzmann constant. In this paper, CeO_{2} is assumed to have a 90% absorption efficiency. Reaction enthalpy of CeO_{2} reduction Q_{reac} in Eq. (4) is around 475 kJ per 1/2 mole of O_{2} [59].
3.2.3 Chemical Reaction
In our study, the mixed chemical kinetic and equilibrium model was applied to simulate the CeO2 reduction of DEM particles. The reaction was set up based on the DEM particle phase and combined the mass transfer with the fluid phase via the former introduced source term ($S_{m,O_2}$).
Based on the thermodynamic equilibrium model, which combines the crystal site model with the standard Gibbs free energy equation [42], the nonstoichiometric coefficient δ is a function of temperature and oxygen partial pressure [31]:
\[ \frac{\delta}{0.35\delta}=\left(\frac{106,000\ Pa}{P_{O_2}}\right)^{0.217}\exp{\left(\frac{195.6{{kJ}/{mol}}}{RT}\right)}\tag{30} \]
The reaction rate in the DEM particle phase is expressed as m_{p} in Eq. (25). To describe the reaction in the particle phase, the CeO_{2} reduction was regarded as the form of particle devolatilization. As the CeO_{2} reduction process is a nonstoichiometric reaction, the user reaction rate method (Eq. (31) is applied instead of the default first order rate method (Eq. (32)) [53].
\[ {\dot{m}}_p=\frac{dm}{dt}=r_{user\ rate}\tag{31} \]
\[ \frac{dm}{dt}=km\tag{32} \]
Based on the study of Ishida et al. [31], the CeO_{2} reduction rate can be expressed as:
\[ \frac{dx}{dt}=k_{red}\left(1x\right)\tag{33} \]
where k_{red} is the reduction rate constant, x is the reaction fraction. The reaction fraction can be conveyed as:
\[ x=\frac{mm_i}{m_fm_i}\tag{34} \]
where m, m_{i}, and m_{f} are timedependent particle mass, initial particle mass, and final particle mass, respectively.
Since Ce is conserved in the CeO_{2} reduction process, the concentration of vacant oxygen ([O_{vac}]) and the concentration of cerium ([Ce]) have the following relation:
\[ \frac{\left[O_{vac}\right]}{\left[Ce\right]}=\delta\tag{35} \]
Based on the above conservation between the concentration of oxygen and the concentration of cerium, the timedependent particle mass m in Eq. (31) can be shown as a function of δ:
\[ m=m_i\frac{\delta n_{CeO_2}M_{O_2}}{2}\tag{36} \]
where $n_{CeO_2}$ is the moles of CeO_{2} and $M_{O_2}$ is the molecular weight of oxygen.
By substituting Eq. (30) and (33) into Eq. (31), the production generation rate is
\[ r_{user\ rate}=\frac{dm}{dt}=\left(m_fm_i\right)\cdot k_{red}\cdot\left(1+\frac{\delta\frac{m_i}{M_{{CeO}_2}}M_{O_2}}{2\left(m_fm_i\right)}\right)\tag{37} \]
The particle devolatilization rate based on Eq. (31) can be expressed as:
\[ {\dot{m}}_p=\left(m_fm_i\right)\cdot k_{red}\cdot\left(1+\frac{\delta\frac{m_i}{M_{{CeO}_2}}M_{O_2}}{2\left(m_fm_i\right)}\right)\tag{38} \]
3.3 Boundary and Initial Conditions
The boundary conditions are summarized in Table 1. A noslip boundary condition is assumed on all fluidsolid interfaces. All the outer wall surfaces are adiabatic. The inlet surfaces are set as velocity inlet condition with the velocity taken from the volume flow rate. The outlets operate under atmospheric pressure. The boundary condition of the partition inner surface is taken as a diffuse radiation flux condition, which is determined by input radiation power. The radiative properties of partition and insulation boundaries depend on the material properties (listed in Table 1).
The entire reduction process can be divided into two sections: 1) preheating and 2) reaction. The initial conditions are distinguished for each section. The preheating process is initiated from room temperature (298 K) and atmospheric pressure (1 atm) with pure Ar purge gas. The particle devolatilization is deactivated. After 20 mins radiation under 1.5 kW power, the preheating process is terminated. In the reaction stage, the initial conditions are taken from the preheating stage, except the inlet oxygen mole fraction is set as 10 ppm.
3.4 Parameters
The material properties of gas components, catalyst, and solid components are listed in Table 2. In this paper, four important parameters are studied to compare the performance via the control variable method. The four parameters can be classified into two groups. One is related to operating conditions, like radiant power and gas flow rate. The other is related to CeO_{2} particle packing, such as void fraction and particle size. To expand the study rationally, a base case is defined with a void fraction of 0.65, a particle size of 5 mm, under radiant power of 3.5 kW in the reaction stage, and the gas flow rate is set to 1.8 L/min (specific parameters listed in Table 3).
Table 3 Studied reactor parameters.
3.5 Numerical Implementation
The transport conservation equations are solved by the finite volume method with approximately 400,000 hexagonal cells in an unstructured mesh using STARCCM+ v. 12.02 [53]. A firstorder implicit unsteady scheme is used for time integration with a time step of 2.0 s. The secondorder segregated flow/energy/species algorithms are applied in the computation of conservation equations. The DEM solver is introduced to solve the governing equations of discrete CeO_{2} particles [60]. The nonstoichiometric CeO_{2} reduction is defined via userdefined field functions as a source term. The transient simulations are performed on a highperformance cluster “FLARE” (6 x (1 Lenovo NeXtScale nx360m5 compute node, 64 GB RAM)).
3.6 Validation
To reduce the influences of mesh number and quality on the simulation results, the mesh independence analysis was applied to this study. Three different cell base sizes were selected to study the average temperature of the catalyst region for 10 minutes preheating process.
The relation between cell numbers and temperatures is shown in Figure 4. With the mesh number increasing, the temperature differences between 400,000 case and 960,000 cases are almost negligible, especially for the gas phase. Therefore, the mesh number reaches a stable level. Considering the accuracy and calculation cost for the transient model, the mesh with 400,000 cells was selected in the following study.
Figure 4 Mesh independence analysis.
The timespatial analysis was applied to this model to ensure the time step is accurate enough to describe the processes in the reactor. Four different time steps (0.5s, 1s, 2s, and 3s) were used in the 20 minutes preheating process. The temperature was monitored by adding a vertical line probe in the center of the catalyst. Figure 5 shows the relation between temperature and position at the moment of 20 minutes under various time steps. The spatialdependent temperature profiles have the same shape and the maximum difference is less than 4%, which is in the tolerable range. Comprehensively considering the transient simulation target, the accuracy, and the calculation cost, 2s time step was chosen in the following simulations.
Figure 5 Timespatial analysis.
The thermodynamic and CeO_{2} reduction models were validated against the experimental results obtained by Bulfin et al. [31], which used 100 W concentrated radiation power on a CeO_{2} pellet (5 mm diameter and 1 mm height) with 65% void fraction. It is visible from Figure 6 that the temperature and nonstoichiometric coefficient of CeO_{2} derived from the numerical model demonstrate comparable values to those derived from the experiment, with slight differences caused by influences of experimental setups and operating limitations [31].
Figure 6 Validation of thermodynamic and chemical reaction models.
4. Modeling Results and Discussions
4.1 Gas Flow Rates
The purge gas is essential in CeO_{2} reduction to remove generated O_{2} from the particle surface to avoid the recombination with reduced CeO_{2}. In this work, the influences of the inlet gas flow rate are discussed. Inlet gas consists of Ar and 10 ppm O_{2}, which was applied to stimulate reduction. Three different flow rates were studied to compare the impacts on the temperatures and the O_{2} evolution rates during the reduction process. The transient temperature profiles of the gas and the particle phase are shown in Figure 7. The results reveal that the temperature differences are too small to distinguish under the given 15 mins reaction time. To clarify the results, temperature distributions between 10 mins and 15 mins are enlarged. From the enlarged temperature profiles, the 1.8 L/min inlet gas flow rate shows relatively higher average temperatures in both phases. As slow gas flow rate takes less cold fluid to the reactive region, the convection losses are smaller than that of large gas flow rate. However, the temperature differences between 1.8 L/min and 3.2 L/min are 5 K and 0.5 K for gas phase and particle phase, respectively. This result reveals that various inlet gas flow rates have a negligible influence on temperature distribution under the current range. Since the input radiant power is large enough to control the temperature distribution, the influence on convection caused by gas flow changes small and may be ignored. The O_{2} evolution rate reveals a similar result (Figure 8). The peak reaction rates between 1.8 L/min and 3.2 L/min is less than 0.01 mL min1 g1 CeO_{2}. Even though the gas flow rate may impact the O_{2} evolution rate, the result is not obvious under current operating conditions.
Figure 7 Temperature profiles of various inlet gas flow rate: (a) gas phase; (b) particle phase.
Figure 8 Oxygen evolution rates of various inlet gas flow rate.
4.2 Solar Radiant Powers
Solar radiant power is a significant factor in CeO_{2} reduction rate since it supplies the required thermal energy to support the endothermic reduction reaction. Four various radiant power levels between 3.0 kW and 4.5 kW were selected to investigate their influence on the reactive region, including incident radiation, temperature distribution, and reduction rate. Figure 9 (a) shows the scattering distribution of the radiant flux of the catalyst region’s inner surface after 15 mins. The Z positive direction is towards the solar receiver window. Higher input power results in higher radiant flux (as high as 4×10^{6} W m^{2}), which also shows that the radiant flux is lower at the ceria top, even though the power is uniformly distributed on the surface. This phenomenon is caused by the diffuse radiation property. In this case, it is assumed that diffuse radiation works in the domains. Since the view factors are evenly distributed in all directions, less radiation accumulates at the top of the ceria domain, which results in a low radiant flux. Additionally, the radiant flux intervals for the different input power levels are nearly equal (consistent with Eq. (15)). The variations of radiant flux on radial direction can be reflected via temperature distribution, as shown in Figure 9 (b). Here, three different levels (top, middle, and bottom) were chosen to investigate the temperature profiles considering both radial and axial influences. At the top level, the temperature is the lowest with a relatively uniform distribution in the radial direction. This is due to the fluid not directly contacting the radiation absorption solid with a lowradiant flux. At the middle and bottom levels, the temperature profile distributions gradually decrease. Since incident radiation decreases due to absorption by solid particles which then radiate radiation, the radiation term in Eq. (18) decreases in the radial direction, which is consistent with a decrease temperature. Due to the high radiant flux to the middle level, the temperature is also higher.
The transient results under different input power are shown in Figure 9 (c) and (d). Figure 9 (c) shows the average temperature variations of the gas and particle phase in their respective domain during the 15 min reaction time. Since particles have higher specific heat and have better radiant properties compared with gases (Ar and O_{2} are weak radiation absorbing species [56]), the particle phase temperatures are higher than the gas phases (difference between gas phase and particle phase is around 300 K). Figure 9 (d) presents that reactive region under a 4.5 kW power input which produces the highest O_{2} formation rate due to its high temperature. Additionally, reaction rates sharply increase after 5 mins, due to temperature exceeds the minimum required reduction temperature (higher than 1000 °C, shown in Figure 9 (c)). After 10 mins, the average gas and particle phase temperatures of the reactive region stabilize, which provides a constant reaction rate. Under 3.5 kW power input, the peak O_{2} evolution rate reaches 0.8 mL min1 g1 CeO_{2}, which is comparable to that reported in the literature [69]. With 4.5 kW power input, the peak O_{2} evolution rate reaches 1.18 mL min1 g1 CeO_{2}. The results show that higher power input provides a higher temperature in the reactive domain and correspondingly leads to higher O_{2} production. However, this observation does not mean the power input can continue to be increased without bound. An optimal power input should also consider solartofuel efficiency. Eq. (7) shows the theoretical solartofuel efficiency ignoring the energy of purge gas [22].
Figure 9 Profiles of various radiant powers: (a) radiant flux distributions along the axial direction; (b) temperature distributions along the radial direction; (c) transient temperature of fluid and particles; (d) oxygen evolution rates.
4.3 CeO_{2} Void Fractions
Besides the influences caused by the operating condition, CeO_{2} morphology is also extremely important in the reduction process. In this work, DEM particles are used to approximate porous media via random injection, which provided similar gas paths among particles. Heat and mass transfer are considered macroscopically, which only occurs on the surface of particles by neglecting particles pores.
In this work, four different void fractions (0.7, 0.65, 0.6, and 0.5) were considered to study temperature distributions and corresponding O_{2} evolution rates. Using the same particle size to fill the reactor’s reactive region, a large void fraction implies that the packing is less dense, which provides less pressure drop in the packed section. Also, a larger void fraction allows incident radiation to penetrate deeper into the reactor radial direction, which results in higher reaction temperature. The resulting temperature distributions for the gas and particle phases are shown in Figure 10 which confirms that tighter (lower void fraction) packing results in lowtemperature distributions for the gas and particle phases consistent with expectation, while higher void fraction (less tightly packed section) results in a higher temperature. Figure 11 shows the timedependent temperature profiles and O_{2} evolution rates. As discussed above, the average temperatures of gas and particle phases are different due to the physical and material properties. In Figure 11 (a), particle temperature increases with an increased void fraction. However, the gas phase temperatures remain essentially constant for void fractions between ε=0.7 and ε=0.65 during the reduction process. From Eq. (18) and (21), the primary heat transfer mechanisms is convection from both the fluidwall exchange and fluidparticle convection. When ε=0.65, the fluidwall convection is greater but has a smaller energy source term for the particle phase. Under the current geometry, the fluidwall convection and fluidparticle source term are balanced, which results in the close temperature distribution with ε=0.70. A similar temperature distribution leads to a similar O_{2} evolution rate in a section with increasing temperature. In Figure 11 (b), the O_{2} evolution rate is almost the same in the initial 7 mins for the void fractions between ε=0.7 and ε=0.65 due to a uniform temperature profile. After that, the temperatures are more or less stable, which reflects a stable oxygen generation rate with a slight difference. Since O_{2} evolution rate is based on per gram CeO_{2}, ε=0.70 has a higher rate caused by the lower mass amount of CeO_{2} under a similar temperature condition. In order words, the catalyst with a void fraction of 0.7 can provide better heat transfer and achieve higher temperature distribution, which leads to a higher O_{2} evolution rate.
Figure 10 Simulated (a) gas phase and (b) particle phase temperature profiles under various CeO_{2} void fractions.
Figure 11 Timedependent variable profiles under different void fractions: (a) average temperature; (b) O_{2} evolution rates.
4.4 CeO_{2} Particle Sizes
CeO_{2} particle size is another important factor that affects the structure of the reactive section. Under the same void fraction, smaller particles provide a more total specific surface area in the reactor. As O_{2} evolution reaction occurs on the surface of CeO_{2} particles, the larger specific surface area is beneficial to the reaction rate. However, smaller particles mean tighter packing, which is adverse to radiative heat transfer. Therefore, the influence of particle sizes may be dual directions. Figure 12 shows the temperature profiles of the gas phase and particle phase with ε=0.65 under various particle diameters. The result reveals that larger particles have advantages of deriving higher temperature since more incident radiation can reach particle surfaces due to the larger interspace between particles, which is shown in Figure 12 (b).
Figure 12 Simulated (a) gas phase and (b) particle phase temperature profiles under different particle sizes.
The transient variations of temperatures and O_{2} evolution rates are presented in Figure 13. From Figure 13 (a), the particle temperature differences between dp=3 mm and dp=6 mm are distinct. The particle diameter dp=3 mm gets a 150 K lower average particle temperature than that of dp=6 mm. Compared with the temperature profile with other cases, dp=3 mm derives the lowest temperature than other cases, which reveals that particle size has a stronger influence on temperature than other factors. Even though the temperature differences between different particle sizes are obvious, the gas phase temperatures seem discrepant. The gas phase temperature seems very close under different conditions, especially for dp=5 mm and dp=6 mm. The result illustrates that the influence of particlegas convection is weaker than that of gaswall convection. In Figure 13 (b), dp=3 mm shows the highest O_{2} evolution rate under these three cases, which is opposite to temperature results. As discussed above, specific surface area and temperature are both significant to the reaction rate. Under the CeO_{2} reduction process, the specific surface area is dominant.
Figure 13 Timedependent variable profiles under different particle sizes: (a) average temperature; (b) O_{2} evolution rates.
5. Conclusions
The DEM method was successfully employed to simulate the porous CeO_{2} structure in the transient CFD modeling of CeO_{2} reduction in a solar thermochemical reactor. Model validation was accomplished according to previously reported experimental data. The effects of the catalyst textures (particle size and void fraction) and process conditions (gas flow rate and radiative power input) on the temperature profiles and reaction rate were investigated. The gas flow rate showed slight influences on both the temperature distribution and reaction rates. The temperature is higher at a lower flow rate leading to a higher reaction rate. The higher solar radiant power input results in higher temperature and O_{2} evolution rate. A larger void fraction of the catalyst is advantageous to improve the thermal performances and CeO_{2} reduction rate. The reaction rate is remarkably enhanced by shrinking the particle size to derive more specific surface area of the catalyst. Further studies will focus on the study of CO_{2} splitting via reduced CeO_{2} and experimental validation of the thermal and reactive performances.
Nomenclature
Symbols 

A 
surface area 
D 
diameter 
$f_b$ 
body force (gravity) 
$F_s$ 
particle surface force (N) 
$F_b$ 
particle body force (N) 
H 
enthalpy (J kg^{1}) 
I 
radiative intensity (W m^{2}) 
I 
identity matrix 
k 
thermal conductivity (W m^{1} K^{1}) 
k_{red} 
reduction rate coefficient 
m_{f} 
final particle mass 
m_{i} 
initial particle mass 
$ \dot{m}_p $ 
rate of mass transfer to particle 
q_{r} 
reradiation (W m^{2}) 
Q 
heat transfer (W) 
Q_{a,p} 
absorption efficiency of particle 
Q_{t} 
surface heat transfer (W) 
s 
distance in Ω direction 
S_{E} 
energy source (W m^{2} or W m^{3}) 
S_{m} 
mass source (mol m^{3} s^{1}) 
t 
time (s) 
T 
temperature (K) 
v 
velocity (m s^{1}) 
x 
reaction fraction 
Y_{i} 
component concentration 
Greek Symbols 

α 
void fraction of fluid phase 
β 
extinction coefficient 
δ 
nonstoichiometric coefficient 
κ_{a} 
absorption coefficient (m^{1}) 
κ_{pa} 
particle absorption coefficient (m^{1}) 
κ_{ps} 
particle scattering coefficient (m^{1}) 
κ_{s} 
scattering coefficient (m^{1}) 
ρ 
density (kg m^{3}) 
σ 
StefanBoltzmann constant 
Ω 
solid angle 
Subscripts 

b 
black body 
conv 
convection 
f 
fluid 
fs 
fluidsolid interphase 
p 
particle 
part 
partition 
pb 
particle black body 
rad 
radiation 
reac 
reaction (reduction) 
s 
insulation 
Abbreviation 

DEM 
Discrete Element Method 
DOM 
Discrete Ordinates Model 
i.d. 
inner diameter 
Nu 
Nusselt number 
o.d. 
outer diameter 
ppm 
part per million 
Pr 
Prandtl number 
Re 
Reynold’s number 
Acknowledgments
We gratefully acknowledge the financial support from the Wayne and Gayle Laufer Foundation.
Author Contributions
Smith, Joseph: Conceptualization, Methodology, Supervision, Reviewing and Editing; Zhang, Han: Software, Data curation, WritingOriginal draft preparation, Visualization, Investigation, Software, Validation, Writing.
Competing Interests
The authors have declared that no competing interests exist.
References
 Energy related CO_{2} emissions [Internet]. Paris: IEA; 2020 [cited date 2020 February 11]. Available from: https://www.iea.org/dataandstatistics/charts/energyrelatedCO_{2}emissions19902019.
 Breyer C, Bogdanov D, Aghahosseini A, Gulagi A, Child M, Oyewo AS, et al. Solar photovoltaics demand for the global energy transition in the power sector. Prog Photovolt Res Appl. 2018; 26: 505523 [CrossRef]
 Zhu T, SwaminathanGopalan K, Stephani K, Ertekin E. Thermoelectric phononglass electroncrystal via ion beam patterning of silicon. Phys Rev B. 2018; 97: 174201. [CrossRef]
 Dincer I, Acar C. A review on clean energy solutions for better sustainability. Int J Energy Res. 2015; 39: 585606. [CrossRef]
 Xu J, Wang F, Lv C, Xie H. Carbon emission reduction and reliable power supply equilibrium based daily scheduling towards hydrothermalwind generation system: A perspective from China. Energy Convers Manag. 2018; 164: 114. [CrossRef]
 Vaughan NE, Gough C, Mander S, Littleton EW, Welfle A, Gernaat DE, et al. Evaluating the use of biomass energy with carbon capture and storage in low emission scenarios. Environ Res Lett. 2018; 13: 044014. [CrossRef]
 Nguyen VN, Blum L. Syngas and synfuels from H_{2}O and CO_{2}: Current status. Chemie Ing Tech. 2015; 87: 354375. [CrossRef]
 Dayton DC, Turk B, Gupta R. Syngas cleanup, conditioning, and utilization. Thermochem Process Biomass Convers into Fuels Chem Power. 2019; 125174. doi: 10.1002/9781119417637.ch5. [CrossRef]
 Bhosale RR. Thermodynamic efficiency analysis of zinc oxide based solar driven thermochemical H_{2}O splitting cycle: Effect of partial pressure of O_{2}, thermal reduction and H_{2}O splitting temperatures. Int J Hydrog Energy. 2018; 43: 1491514924. [CrossRef]
 Bachirou GL, Shuai Y, Zhang J, Huang X, Yuan Y, Tan H. Syngas production by simultaneous splitting of H_{2}O and CO_{2} via iron oxide (Fe3O4) redox reactions under highpressure. Int J Hydrog Energy. 2016; 41: 1993619946. [CrossRef]
 Wu S, Zhou C, Doroodchi E, Nellore R, Moghtaderi B. A review on hightemperature thermochemical energy storage based on metal oxides redox cycle. Energy Convers Manag. 2018; 168: 421453. [CrossRef]
 Scheffe JR, Steinfeld A. Oxygen exchange materials for solar thermochemical splitting of H_{2}O and CO_{2}: A review. Mater Today. 2014; 17: 341348. [CrossRef]
 Chueh WC, Falter C, Abbott M, Scipio D, Furler P, Haile SM, et al. Highflux solardriven thermochemical dissociation of CO_{2} and H_{2}O using nonstoichiometric ceria. Science. 2010; 330: 17971801. [CrossRef]
 Romero M, Steinfeld A. Concentrating solar thermal power and thermochemical fuels. Energy Environ Sci. 2012; 5: 92349245. [CrossRef]
 Muhich CL, Ehrhart BD, Al‐Shankiti I, Ward BJ, Musgrave CB, Weimer AW. A review and perspective of efficient hydrogen generation via solar thermal water splitting. Wiley Interdiscip Rev Energy Environ. 2016; 5: 261287. [CrossRef]
 Charvin P, Abanades S, Flamant G, Lemort F. Twostep water splitting thermochemical cycle based on iron oxide redox pair for solar hydrogen production. Energy. 2007; 32: 11241133. [CrossRef]
 Dizaji HB, Hosseini H. A review of material screening in pure and mixedmetal oxide thermochemical energy storage (TCES) systems for concentrated solar power (CSP) applications. Renew Sust Energy Rev. 2018; 98: 926. [CrossRef]
 Bulfin B, Vieten J, Agrafiotis C, Roeb M, Sattler C. Applications and limitations of two step metal oxide thermochemical redox cycles: A review. J Mater Chem A. 2017; 5: 1895118966. [CrossRef]
 Otsuka K, Hatano M, Morikawa A. Hydrogen from water by reduced cerium oxide. J Catal. 1983; 79: 493496. [CrossRef]
 Abanades S, Flamant G. Thermochemical hydrogen production from a twostep solardriven watersplitting cycle based on cerium oxides. Sol Energy. 2006; 80: 16111623. [CrossRef]
 Petrasch J, Klausner J. Integrated solar thermochemical cycles for energy storage and fuel production. Wiley Interdiscip Rev Energy Environ. 2012; 1: 347361. [CrossRef]
 Bader R, Chandran RB, Venstrom LJ, Sedler SJ, Krenzke PT, De Smith RM. Design of a solar reactor to split CO_{2} via isothermal redox cycling of ceria. J Sol Energy Eng. 2015; 137. [CrossRef]
 Furler P, Scheffe JR, Steinfeld A. Syngas production by simultaneous splitting of H_{2}O and CO_{2} via ceria redox reactions in a hightemperature solar reactor. Energy Environ Sci. 2012; 5: 60986103. [CrossRef]
 Furler P, Scheffe J, Gorbar M, Moes L, Vogt U, Steinfeld A. Solar thermochemical CO_{2} splitting utilizing a reticulated porous ceria redox system. Energy Fuels. 2012; 26: 70517059. [CrossRef]
 Venstrom LJ, Petkovich N, Rudisill S, Stein A, Davidson JH. The effects of morphology on the oxidation of ceria by water and carbon dioxide. J Sol Energy Eng. 2012; 134: 11005. [CrossRef]
 Zoller S, Koepf E, Roos P, Steinfeld A. Heat transfer model of a 50kW solar receiverreactor for thermochemical redox cycling using cerium dioxide. J Sol Energy Eng. 2019; 141: 21014. [CrossRef]
 Sedighi M, Padilla RV, Taylor RA, Lake M, Izadgoshasb I, Rose A. Hightemperature, pointfocus, pressurised gasphase solar receivers: A comprehensive review. Energy Convers Manag. 2019; 185: 678717. [CrossRef]
 de la Calle A, Bayon A. Annual performance of a thermochemical solar syngas production plant based on nonstoichiometric CeO_{2}. Int J Hydrog Energy. 2019; 44: 14091424. [CrossRef]
 Groehn AJ, Lewandowski A, Yang R, Weimer AW. Hybrid radiation modeling for multiphase solarthermal reactor systems operated at hightemperature. Sol Energy. 2016; 140: 130140. [CrossRef]
 Bulfin B, Call F, Lange M, Lubben O, Sattler C, PitzPaal R, et al. Thermodynamics of CeO_{2} thermochemical fuel production. Energy Fuels. 2015; 29: 10011009. [CrossRef]
 Bulfin B, Lowe AJ, Keogh KA, Murphy BE., Lubben O, Krasnikov SA, et al. Analytical model of CeO_{2} oxidation and reduction. J Phys Chem C. 2013; 117: 2412924137. [CrossRef]
 Zinkevich M, Djurovic D, Aldinger F. Thermodynamic modelling of the cerium–oxygen system. Solid State Ion. 2006; 177: 9891001. [CrossRef]
 Ishida T, Gokon N, Hatamachi T, Kodama T. Kinetics of thermal reduction step of thermochemical twostep water splitting using CeO_{2} particles: Masterplot method for analyzing nonisothermal experiments. Energy Procedia. 2014; 49: 19701979. [CrossRef]
 ValadesPelayo PJ, VillafánVidales HI, RomeroParedes H, ArancibiaBulnes CA. Modeling of a tubular solar reactor for continuous reduction of CeO_{2}: The effect of particle size and loading on radiative heat transfer and conversion. Chem Eng Sci. 2017; 162: 7787. [CrossRef]
 ValleHernández J, RomeroParedes H, PachecoReyes A. Modeling of the steam hydrolysis in a twostep process for hydrogen production by solar concentrated energy. AIP Conf Proc. 2017; 1850: 100016. [CrossRef]
 ValleHernández J, RomeroParedes H, ArancibiaBulnes CA, VillafanVidales HI, EspinosaParedes G. Modeling of a CeO_{2} thermochemistry reduction process for hydrogen production by solar concentrated energy. AIP Conf Proc. 2016; 1734: 120008. [CrossRef]
 Arifin D, Weimer AW. Kinetics and mechanism of solarthermochemical H_{2} and CO production by oxidation of reduced CeO_{2}. Sol Energy. 2018; 160: 178185. [CrossRef]
 McDaniel AH, Ambrosini A, Coker EN, Miller JE, Chueh WC, O’Hayre R, et al. Nonstoichiometric perovskite oxides for solar thermochemical H_{2} and CO production. Energy Procedia. 2014; 49: 20092018. [CrossRef]
 Scheffe JR, McDaniel AH, Allendorf MD, Weimer AW. Kinetics and mechanism of solarthermochemical H_{2} production by oxidation of a cobalt ferritezirconia composite. Energy Environ Sci. 2013; 6: 963973. [CrossRef]
 Kröger FA, Vink HJ. Relations between the concentrations of imperfections in crystalline solids. Solid State Phys. 1956; 3: 307435. [CrossRef]
 Scheffe JR, Steinfeld A. Thermodynamic analysis of ceriumbased oxides for solar thermochemical fuel production. Energy Fuels. 2012; 26: 19281936. [CrossRef]
 Bader R, Gampp L, Breuille T, Haussener S, Steinfeld A, Lipiński W. Unsteady Radiative heat transfer model of a ceria particle suspension undergoing solar thermochemical reduction. J Thermophys Heat Transf. 2019; 33: 6377. [CrossRef]
 Haeussler A, Abanades S, Julbe A, Jouannaux J, Cartoixa B. Solar thermochemical fuel production from H_{2}O and CO_{2} splitting via twostep redox cycling of reticulated porous ceria structures integrated in a monolithic cavitytype reactor. Energy. 2020; 117649. [CrossRef]
 Haeussler A, Abanades S, Julbe A, Jouannaux J, Drobek M, Ayral A, et al. Remarkable performance of microstructured ceria foams for thermochemical splitting of H_{2}O and CO_{2} in a novel hightemperature solar reactor. Chem Eng Res Des. 2020; 156: 311323. [CrossRef]
 Parthasarathy P, Le Clercq P. Heat transfer simulation in a high temperature solar reactor. Energy Procedia. 2015; 69: 18101818. [CrossRef]
 Kyrimis S, Le Clercq P, Brendelberger S. 3D modelling of a solar thermochemical reactor for MW scalingup studies. AIP Conf Proc. 2019; 2126: 180013. [CrossRef]
 Ma H, Xu L, Zhao Y. CFDDEM simulation of fluidization of rodlike particles in a fluidized bed. Powder Technol. 2017; 314: 355366. [CrossRef]
 Zhang G, Gutierrez M, Li M. A coupled CFDDEM approach to model particlefluid mixture transport between two parallel plates to improve understanding of proppant micromechanics in hydraulic fractures. Powder Technol. 2017; 308: 235248. [CrossRef]
 Bellan S, Matsubara K, Cho HS, Gokon N, Kodama T. A CFDDEM study of hydrodynamics with heat transfer in a gassolid fluidized bed reactor for solar thermal applications. Int J Heat Mass Transf. 2018, 116: 377392. [CrossRef]
 Bellan S, Kodama T, Matsubara K, Gokon N, Cho HS, Inoue K. Thermal performance of a 30kW fluidized bed reactor for solar gasification: A CFDDEM study. Chem Eng J. 2019; 360: 12871300. [CrossRef]
 Morris AB, Ma Z, Pannala S, Hrenya CM. Simulations of heat transfer to solid particles flowing through an array of heated tubes. Sol Energy. 2016, 130: 101115. [CrossRef]
 STARCCM+ Theory Guide v. 12.06. Siemens AG. 2018.
 STARCCM+ User Manual v. 12.06. 2018.
 Furler P, Steinfeld A. Heat transfer and fluid flow analysis of a 4kW solar thermochemical reactor for ceria redox cycling. Chem Eng Sci. 2015; 137: 373383. [CrossRef]

Khalsa AK, De Andrade D. Thermal reflectivity in solar collector prototype. Senior Design Project; 2008. Available from: http://solar.sdsu.edu/Senior%20Projects/Solar_Trough/Solar_Trough_Design
_Khalsa_DeAndrade_Final_Report.pdf.  Brewster MQ. Thermal radiative transfer and properties. Hoboken: John Wiley Sons; 1992.
 Hanley HJ, McCarty RD, Sengers JV. Viscosity and thermal conductivity coefficients of gaseous and liquid oxygen. 1974. Available from: https://core.ac.uk/download/pdf/42893032.pdf.
 Ranz WE, Marshall WR. Evaporation from drops. Chem Eng Prog. 1952; 48: 141146.
 Marxer D, Furler P, Takacs M, Steinfeld A. Solar thermochemical splitting of CO_{2} into separate streams of CO and O_{2} with high selectivity, stability, conversion, and efficiency. Energy Environ Sci. 2017; 10: 11421149. [CrossRef]
 Cundall PA, Strack OD. A discrete numerical model for granular assemblies. Geotechnique. 1979; 29: 4765. [CrossRef]
 Touloukian YS. Thermophysical properties of high temperature solid materials. Volume 4. Oxides and their solutions and mixtures. Part 1. Simple oxygen compounds and their mixtures, thermophysical and electronic properties information analysis center lafayette in. 1966.
 Yaws CL. Thermophysical properties of chemicals and hydrocarbons. Norwich, NY: William Andrew; 2008.
 STARCCM+ Material Databases. Siemens AG. 2017.
 Richet P, Bottinga Y, Denielou L, Petitet JP, Tequi C. Thermodynamic properties of quartz, cristobalite and amorphous SiO_{2}: Drop calorimetry measurements between 1000 and 1800 K and a review from 0 to 2000 K. Geochim Cosmochim Acta. 1982; 46: 26392658. [CrossRef]
 Andre S, Degiovanni A. A theoretical study of the transient coupled conduction and radiation heat transfer in glass: Phonic diffusivity measurements by the flash technique. Int J Heat Mass Transf. 1995; 38: 34013412. [CrossRef]
 Nicolau VP, Borges MT, Silva LM. Radiative properties of glass and coatings. 3rd Eur Conf Energy Perform Indoor Clim Build. 2002; 835840.
 McBride BJ. Coefficients for calculating thermodynamic and transport properties of individual species. Washington: NASA Langley Research Center; 1993.
 Zimmermann D. Flow modeling of a solar thermogravimeter. Master Thesis ETH Zurich. 2012.
 Venstrom LJ, De Smith RM, Chandran RB, Boman DB, Krenzke PT, Davidson JH. Applicability of an equilibrium model to predict the conversion of CO_{2} to CO via the reduction and oxidation of a fixed bed of cerium dioxide. Energy Fuels. 2015; 29: 81688177. [CrossRef]