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 Alonso-Montesinos

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.

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₂ 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₂ emission [2,3,4,5,6]. Since CO₂ 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₂ and H₂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 Fischer-Tropsch technology [8]. Solar energy is essentially a green technology with inexhaustible power production potential. Under the concentrated energy from the sun, CO₂ and H₂O can be converted to syngas using proper metal oxide redox pairs such as ZnO/Zn [9], Fe₃O₄/FeO [10], and many more advanced couples [11]. The appearance of producing syngas via solar energy provides a brand-new idea for both solar energy application and chemical fuel production.

The two-step solar thermochemical process is a redox cycle based on metal oxide redox. This chemical process is represented as the high-temperature reduction (see Eq. (1)) and low-temperature oxidation (see Eq. (2) and (3)).

High-temperature reduction:

\[ M_xO_y\ \rightarrow\ \ M_xO_{y-\delta}+\ \frac{\delta}{2}\ O_2\tag{1} \]

Low-temperature 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₂/H₂O react with metal or reduced metal oxide and lose oxygen to form CO/H₂. 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 two-step 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 two-step solar thermochemical redox process can be divided into two groups, volatile and non-volatile. 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₂/SnO, GeO₂/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]. Non-volatile 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 non-volatile 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 non-stoichiometric CeO₂ redox pair for solar thermochemical reactions flourished in the last decade. In 2006, Flamant et al. [20] proposed the application of CeO₂/Ce₂O₃ redox pair to split water. The reduction of CeO₂ was proceeded at 2000 °C, 100~200 mbar, followed by the oxidation step in the temperature range of 400-600 °C. Since reduction temperature is high enough to transfer CeO₂ to a molten state, which strongly reduces the catalyst activity and quality, the study of CeO₂ was focused on how to avoid the phase transfer under a high reduction temperature. Non-stoichiometric CeO₂ (CeO_2-δ) was proposed to solve the issue of melting at high temperature. Based on the two-step 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₂ splitting ($ CO_2\to CO+\frac{1}{2}O_2 $) and H₂O splitting ($ H_2 O\to H_2+\frac{1}{2}O_2 $) are 283 kJ per mole of CO₂ and 242 kJ per mole of H₂O, respectively [21,22]. At 1500 and 0.1 mbar O₂ partial pressure, the enthalpy change of ceria reduction is around 475 kJ per half mole of O₂ [21]. Therefore, the reaction enthalpies to generate one mole of CO and H₂ are estimated -192 kJ/mol and -233 kJ/mol.

The solar-to-fuel efficiency (see Eq. (7)) is used to show the reaction performance [22].

\[ \eta_{solar-to-fuel}=\frac{\eta_{fuel}HHV}{Q_{solar}}\tag{7} \]

Chueh et al. [13] presented the application of porous monolithic ceria on the two-step solar thermochemical CO₂/H₂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 solar-to-fuel efficiency. Furler et al. [23] proposed using porous ceria felt as the catalyst to support solar CO₂/H₂O splitting. The O₂ and fuel production are 2.89±0.27 mL g-1 CeO₂ and 5.88±0.43 mL g-1 CeO₂. Furler et al. [24] also proposed a reticulated porous ceria (RPC) foam structure for redox reactions, which derived a 3.53% peak solar-to-fuel efficiency. The efficiency improved more than four times compared with former porous monolithic ceria. Venstrom et al. [25] presented a novel morphology of CeO₂, three-dimensional-ordered macroporous (3DOM) CeO₂, 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 non-ordered mesoporous ceria. Oliveira et al. [26] further studied the 3DOM ceria structure for solar thermochemical CO₂ splitting in a tubular reactor. The result showed that the maximum fuel production rate is three times higher than the dual-scale porosity ceria reticulated porous foam.

Based on the difference of heat transfer type, solar thermochemical reactors can be classified as indirectly-irradiated reactors and directly-irradiated reactors. For indirectly-irradiated 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 directly-irradiated 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 directly-irradiated reactors, the reactor can achieve a higher temperature. To fully utilize the solar radiation, the cavity-receiver reactor is commonly used for directly-irradiated solar thermochemical processes. In a cavity-receiver reactor, the space provided in the cavity provides or enhances radiation to the catalyst inner surface, which forms an approximate blackbody. Different types of cavity-receiver reactors have been discussed and used over the past several years to prove that direct irradiation has a huge advantage on the two-step solar H₂O/CO₂ splitting process [27].

To better predict and validate the CeO₂ supported solar thermochemical two-step 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₂ two-step redox reactions. Calle et al. [28], Groehn et al. [29], Bulfin et al. [30] introduced the relation between non-stoichiometric 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 receiver-reactor, which obtained the solar-to-fuel efficiency larger than 10% with dual-scale 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 solid-state 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 solid-state 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₂ 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₂ reduction performances in a particle suspension system. Haeussler et al. [43,44] studied the two-step redox cycling of reticulated porous ceria structures in a monolithic cavity-type 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 solar-to-fuel 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 scale-up solar thermochemical reactor geometry. The results showed that the scale-up geometry has a benefit for the solar-to-fuel efficiency at the beginning of the reduction but limited due to the thickness of CeO₂ RPC.

Gas-solid performance plays a vital role in the numerical simulation of solar thermochemical processes. The Eulerian-Lagrangian (CFD-DEM) 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 high-temperature 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 CFD-DEM method and the results showed high agreement between the experimental and simulation results [49,50]. Morris et al. used CFD-DEM 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 CFD-DEM method, very few works consider the reactions in the solar thermochemical two-step 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 three-dimensional (3D) transient multi-physics computational fluid dynamics model is applied to study heat and mass transfer of a high-temperature CeO₂ 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₂ 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₂ reduction in the novel partition-cavity 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 stainless-steel 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) cross-section.

3. Approach

3.1 Reactor Model

Figure 3 shows the solar reactor modeling geometry in top view and cross-section. Since the reactor has a symmetric structure, only one-eighth 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 stainless-steel 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₂) 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₂O₃-SiO₂ 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₂ is assumed to be a black body radiator and Al₂O₃-SiO₂ 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 non-reactive regions and reactive regions. Non-reactive 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₂O₃-SiO₂ 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₂ particles. To set up the link between the particles and fluid phases, a two-way coupling method is introduced that solves the heat and mass transfer conservation equations in both phases [52].

3.2.1 Non-reactive 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 re-radiation term. The input radiation and the re-radiation 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₂ are low radiative materials, the radiative properties κ_a and κ_s are assumed to equal zero [56]. Re-radiation 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 core-reactive region in the reactor consisting of two phases: gas mixture fluid phase and CeO₂ particle Lagrangian phase. The fluid continuous phase should consider the interaction between gas mixture and CeO₂ 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_s-T_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₂ 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(T-T_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 Ranz-Marshall 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(G-4\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 Stefan-Boltzmann constant. In this paper, CeO₂ is assumed to have a 90% absorption efficiency. Reaction enthalpy of CeO₂ reduction Q_reac in Eq. (4) is around 475 kJ per 1/2 mole of O₂ [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 non-stoichiometric 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₂ reduction was regarded as the form of particle devolatilization. As the CeO₂ reduction process is a non-stoichiometric 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₂ reduction rate can be expressed as:

\[ \frac{dx}{dt}=k_{red}\left(1-x\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{m-m_i}{m_f-m_i}\tag{34} \]

where m, m_i, and m_f are time-dependent particle mass, initial particle mass, and final particle mass, respectively.

Since Ce is conserved in the CeO₂ 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 time-dependent 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₂ 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_f-m_i\right)\cdot k_{red}\cdot\left(1+\frac{\delta\frac{m_i}{M_{{CeO}_2}}M_{O_2}}{2\left(m_f-m_i\right)}\right)\tag{37} \]

The particle devolatilization rate based on Eq. (31) can be expressed as:

\[ {\dot{m}}_p=-\left(m_f-m_i\right)\cdot k_{red}\cdot\left(1+\frac{\delta\frac{m_i}{M_{{CeO}_2}}M_{O_2}}{2\left(m_f-m_i\right)}\right)\tag{38} \]

3.3 Boundary and Initial Conditions

The boundary conditions are summarized in Table 1. A no-slip boundary condition is assumed on all fluid-solid 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).

Table 1 Boundary conditions.

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₂ 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 2 Material properties.

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 STAR-CCM+ v. 12.02 [53]. A first-order implicit unsteady scheme is used for time integration with a time step of 2.0 s. The second-order 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₂ particles [60]. The non-stoichiometric CeO₂ reduction is defined via user-defined field functions as a source term. The transient simulations are performed on a high-performance 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 pre-heating 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 time-spatial 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 pre-heating 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 spatial-dependent 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 Time-spatial analysis.

The thermodynamic and CeO₂ reduction models were validated against the experimental results obtained by Bulfin et al. [31], which used 100 W concentrated radiation power on a CeO₂ pellet (5 mm diameter and 1 mm height) with 65% void fraction. It is visible from Figure 6 that the temperature and non-stoichiometric coefficient of CeO₂ 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₂ reduction to remove generated O₂ from the particle surface to avoid the recombination with reduced CeO₂. In this work, the influences of the inlet gas flow rate are discussed. Inlet gas consists of Ar and 10 ppm O₂, which was applied to stimulate reduction. Three different flow rates were studied to compare the impacts on the temperatures and the O₂ 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₂ 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 min-1 g-1 CeO₂. Even though the gas flow rate may impact the O₂ 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₂ 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⁶ 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 low-radiant 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₂ 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₂ 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₂ evolution rate reaches 0.8 mL min-1 g-1 CeO₂, which is comparable to that reported in the literature [69]. With 4.5 kW power input, the peak O₂ evolution rate reaches 1.18 mL min-1 g-1 CeO₂. The results show that higher power input provides a higher temperature in the reactive domain and correspondingly leads to higher O₂ production. However, this observation does not mean the power input can continue to be increased without bound. An optimal power input should also consider solar-to-fuel efficiency. Eq. (7) shows the theoretical solar-to-fuel 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₂ Void Fractions

Besides the influences caused by the operating condition, CeO₂ 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₂ 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 low-temperature 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 time-dependent temperature profiles and O₂ 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 fluid-wall exchange and fluid-particle convection. When ε=0.65, the fluid-wall convection is greater but has a smaller energy source term for the particle phase. Under the current geometry, the fluid-wall convection and fluid-particle source term are balanced, which results in the close temperature distribution with ε=0.70. A similar temperature distribution leads to a similar O₂ evolution rate in a section with increasing temperature. In Figure 11 (b), the O₂ 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₂ evolution rate is based on per gram CeO₂, ε=0.70 has a higher rate caused by the lower mass amount of CeO₂ 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₂ evolution rate.

Figure 10 Simulated (a) gas phase and (b) particle phase temperature profiles under various CeO₂ void fractions.

Figure 11 Time-dependent variable profiles under different void fractions: (a) average temperature; (b) O₂ evolution rates.

4.4 CeO₂ Particle Sizes

CeO₂ 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₂ evolution reaction occurs on the surface of CeO₂ 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₂ 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 particle-gas convection is weaker than that of gas-wall convection. In Figure 13 (b), dp=3 mm shows the highest O₂ 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₂ reduction process, the specific surface area is dominant.

Figure 13 Time-dependent variable profiles under different particle sizes: (a) average temperature; (b) O₂ evolution rates.

5. Conclusions

The DEM method was successfully employed to simulate the porous CeO₂ structure in the transient CFD modeling of CeO₂ 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₂ evolution rate. A larger void fraction of the catalyst is advantageous to improve the thermal performances and CeO₂ 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₂ splitting via reduced CeO₂ and experimental validation of the thermal and reactive performances.

Nomenclature

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, Writing-Original draft preparation, Visualization, Investigation, Software, Validation, Writing.

Competing Interests

The authors have declared that no competing interests exist.