Optimization of Inert Gas Feeding Strategy in a Fixed-Bed Reactor for Efficient Water Splitting Via Solar-Driven Thermal Reduction of Nonstoichiometric CeO₂

Hydrogen and carbon monoxide are clean and carbon-neutral alternative fuels that can be produced from renewable energy sources such as solar energy. There are three main paths for producing H₂/CO from H₂O/CO₂ using solar energy: the solar electrochemical, photochemical, and solar thermochemical routes [1]. In the solar thermochemical route, solar radiation of the entire solar spectrum is concentrated to obtain the high temperature required to drive the thermochemical cycles [2]. Through thermochemical cycles, the products of water splitting, H₂/CO and O₂, are generated in different steps of the thermochemical cycles, bypassing the fuel/oxygen separation issue. Meanwhile, the thermochemical cycles are operated at temperatures significantly lower than thermolysis (>2500 K) [3]. Thermochemical cycles with different numbers of reaction steps are available for water splitting [4–6]. The two-step redox cycles are most promising due to their simplicity and high theoretical efficiency [1]. Such systems employ concentrated solar radiation as a high-temperature heat source to drive an endothermic reduction reaction [7]. The reduced metal-oxide is then oxidized with the water/CO₂ and syngas is produced [8,9].

A wide variety of metal-oxide redox pairs were investigated to conduct the two-step redox cycle, such as Fe₃O₄/FeO, TiO₂/TiO_x, Mn₃O₄/MnO, Co₃O₄/CoO, and ZnO/Zn [10–15]. The CeO₂ (ceria) offers a major promise for efficient fuel generation due to its fast kinetics and crystallographic stability [6,16]. Various experimental studies have investigated the thermodynamics of reactive cerium-based oxides [17]. High-temperature thermodynamic data of undoped and doped ceria have become available for analyzing the thermodynamic processes involving ceria [8,17–21].

Many solar reactor concepts have been proposed for metal-oxide redox cycles including integrated receiver reactors with moving or stationary fixed beds of redox materials, separated reactors with fluidized beds, and others, e.g., aerosol flow reactors [22–32]. The majority of the reactor systems demonstrated so far were based on the fixed-bed type due to the design simplicity and convenient operation [33]. Three generations of solar cavity receiver/reactors containing porous monolithic ceria with multi-scale porosity were developed, achieving consecutive redox cycles for solar fuel production at different power scales [16,22,34–36]. A solar-to-fuel energy efficiency of 5.25% was reported in 2017, which is the current world record for the thermochemical splitting of CO₂ into CO and O₂ [37].

Solar-driven reduction is the crucial phase of thermochemical H₂/CO production, influencing the overall thermal efficiency. According to experimental data, thermodynamically, ceria reduction is favored at high temperatures and low oxygen partial pressure [18]. Prior studies have demonstrated that cavity-shape receivers enable to reach a high temperature as exposed in high-flux solar irradiation (∼3000 suns). However, as ceria undergoes a phase change at temperatures above 2100 K, the maximum reduction temperature is limited [38,39]. On the other hand, two methods are commonly used to achieve a low oxygen partial pressure: vacuum pumping and inert gas sweeping. Here we focused on the latter one because of its advantages in solid-phase heat recovery. Inert gases, e.g., nitrogen or argon, were used as oxygen scavengers to sweep the reactant bodies and carry away released oxygen during the reduction process. References [16,34,35] showed the experimental steps for sweeping strategies of inert gases. The solar reduction is a complex mass and heat transfer process with intense nonlinear temperature changes due to reduction rates, cavity heat emissions, heat losses by gas flows, etc. Optimal solutions with regard to the nonlinear process have not been explored sufficiently yet. First, the absence of a general definition of the characteristic parameter for optimal reaction duration impeded understanding the mechanism for this kind of solar-driven fixed-bed thermal reduction processes. Optimal duration should be determined by weighing the increase of chemical enthalpies and heat losses over time. Also, the global optimization issue related to variable inert gas feeding strategies was lack of comprehensive analysis to maximize overall thermal efficiency. These research gaps are the main focus of this paper. As the genetic algorithm (GA) constitutes a broadly applicable base for stochastic and global optimization [40] and multi-variate optimizing problems [41], we employed it here to address the optimization issue of N₂ feeding strategy. The novel N₂ feeding strategy presented in this study outperforms the state-of-the-art strategies and results in an unprecedented thermal efficiency.

The paper presents the transient model for a cavity fixed-bed reduction system using a lumped parameter method based on energy balance equations in Sec. 2. Section 3 presents (1) the results of the reduction heat loss and thermal efficiency and the best reducing time under different scales of constant flowrates of sweeping gas and (2) the optimal solutions of variable flowrates found by the genetic algorithm under specific conditions. A conclusion is given in Sec. 4.

The schematics of the initial and final phases of the solar-driven ceria-based fixed-bed reduction system are illustrated in Fig. 1. The system consists of a cavity reactor and gas-phase heat exchangers (HXers). For the initial phase, transient thermodynamic states in the system are denoted by numbered circles ①–④ as SP1–SP4. The thin solid arrows and thick hollowed arrows indicate the flow of mass and energy, respectively. The rectangular dash border denotes the system boundary. Subscripts of temperatures represent the phase and state points (SPs), e.g., T_0,2 denotes the temperature at the initial phase point and SP2. Purified N₂ with low oxygen partial pressure $(PO2,0=10−5)$ is used. It enters the system at ambient temperature and pressure (SP1) in a transient flowrate $m˙N2$⁠, preheated by the HXer to SP2 before entering the cavity. A certain amount of nonstoichiometric cerium oxide is filled in the cavity. It is reduced with concentrated solar irradiation and sweeping gas. The released O₂ is carried away by the N₂ flow, which is reheated to SP3. After passing through the HXer again, the mixture gases leave the system at SP4. The following assumptions and baseline parameters are used in the system model:

1. The geometries of the cavity reactor, HXer, and pipes are not considered.

2. Direct normal irradiation is set to 1000 W/m² according to the solar radiation resource available in the western region of China [42].

3. Solar flux concentration ratio is set to 3000² unless otherwise specified [3,43].

4. Heat losses from conduction are lumped into the Q_cond with a loss fraction F_C set to 0.2 [38,44].

5. Effectiveness of gas-phase heat recovery ɛ refers to 95% according to current manufacturing levels of heat exchangers [44–46].

6. Reduction process is in a quasi-equilibrium thermodynamic state.

7. Ceria is treated as homogeneous and isotropic and gases are ideal.

8. Modeling of oxidation is assumed as complete water splitting at constant T_ox equal to 1100 K.

In our case, ceria is reduced from the initial phase of nonstoichiometry (δ_int) to the final one (δ_fin) over a certain period. For any moment, δ equilibrates at corresponding T and $PO2$⁠, obtained by a certain number of iterations with the unconstrained nonlinear optimization method. Equations (11) and (12) were discretized using the backward Euler method and incorporated it in the matlab program. The initial temperature (T_int) and oxygen partial pressure $(PO2,int)$ are set to 1100 K and 10⁻². Since the time scale of reduction varies with the mass of ceria, the specific time-step is defined as Δκ = Δt/m_ceria, set to 1 s/kg_ceria based on the trade-off between accuracy and CPU time needed. The baseline model parameters and expressions are given in Table 1.

First, the conventional parametric analysis method (PAM) is used to search the optimal solutions for the highest thermal efficiencies (η_th) under the condition of constant N₂ flowrates $(m˙N2,const)$⁠. The term $m˙N2,const$ varies in the range of 4 × 10⁻⁵ to 0.0067 kg/s. Note that the reduction duration (τ) is proportional to the mass of ceria (m_ceria) under the same conditions. For simplicity, the specific duration (κ = τ/m_ceria) is used instead of τ in the parameter analysis. The optimized specific duration (κ_pot), corresponding to the maximum thermal efficiency (η_th,max), is found when ∂η_th/∂κ = 0, for given $m˙N2,const$⁠.

The genetic algorithm (GA) approach is used to seek the global optimization solution regarding the feeding strategy of variable inert gas flows. The GA is a probabilistic search procedure focusing on large spaces involving states that can be represented by strings based on the mechanics of natural selection and natural genetics [48]. It consists of four procedures: initialization, selection, reproduction, and termination. GA starts with the initial population, a set of random solutions represented by chromosomes. At each generation, each chromosome is evaluated to be selected or not according to the fitness function composed of one or multiple objectives. Then survivors randomly mate and produce the next generation of chromosomes while crossover and mutation randomly occur. Because chromosomes with high fitness values have a high probability of survival, the average fitness values of the new generation are higher than the old generation. The process is repeated until the stopping criteria, e.g., the limit of generation numbers or the function tolerance, are satisfied.

Note that Eq. (16) is derived from Eqs. (11) and (14) to a limiting case, i.e., T = T_i, dT/dt = 0, and dδ/dt = 0, when t = ∞. It indicates a maximum x_i allowed for the ith stage, ensuring that the target temperature (T_i) is attainable at the end of the ith stage.

Population sizes are set to 50 for n ≤ 5 and to 200 for n > 5. The creation function of initial population is specified as constraint dependent according to Eq. (16). Raw fitness scores returned by the fitness function are used to rank all individuals as input values for the stochastic uniform selection function. The number of individuals guaranteed to survive to the next generation is specified to 5% of population sizes, and the crossover fraction is set to 0.8. Mutation and crossover functions are specified as adaptive feasible and scattered. Optimizations are terminated when the average changes in the fitness values are less than 10⁻⁶.

Prior to the optimization, the fitness function needs to be validated. Numerical and experimental data of the reducing processes via a fixed-bed of ceria are quantitatively compared for the validation [50]. First, the predicted results of O₂ production rates $(rO2)$⁠, by Venstrom’s and the model developed in this study, are compared under the same operating conditions. The ceria is heated from 1013 K to 1748 K at a rate of 100 K/min followed by a steady state for 4 min under the same temperature and a constant N₂ sweeping gas flow of 900 mL/(min · g_ceria). The mass of ceria is 1.0245 g. The initial δ is assumed as 10⁻⁵. Figure 3 shows that the numerical and experimental data are in good agreement. In addition, a horizontal reactor by an infrared imaging furnace was setup for achieving an isothermal condition for the reduction reaction reported in Venstrom’s work [50]. The reactor consists of a packed bed of ceria porous particles. Flowrates of the generated O₂ were measured in real time using a Raman laser gas analyzer. Figures 4(a) and 4(b) compare the numerical results of the O₂ production rate during reduction to the experimental data corresponding to N₂ gas flowrates of 150 ml/(min g) and 300 ml/(min g) at constant reduction temperature of 1773 K. The shaded regions in the plots indicate the range of rates predicted by the model as a result of uncertainty in the initial δ (0.021 ± 0.002). In general, a good agreement in the numerical and experimental data indicates that our model can accurately predict the transient thermal reducing processes of fixed-bed ceria.

Figure 5 shows the bed temperature (T), the oxygen partial pressure $(PO2)$⁠, the thermal efficiency (η_th), and the nonstoichiometric number of ceria (δ) as a function of specific duration (κ) under the conditions of a constant N₂ mass flow $(m˙N2,const=0.005kg/s)$ and different concentration ratios (C varying from 1500 to 3000). As shown in Fig. 5(a), the bed temperature increases monotonously with time for a given solar concentration ratio. Temperature rises from 1100 K to 1816 K in 2 × 10⁴ s/kg under the condition of $m˙N2,const=0.005kg/s$ and C = 1500. With a greater C, the final temperature (T_fin) rises due to a lower cavity radiative loss. Temperatures 1907 K, 1978 K, and 2036 K are obtained at the end of reduction, corresponding to C = 2000, 2500, and 3000, respectively. The oxygen partial pressure increases first and then decreases with time. The peak value of $PO2$ increases with improved C. That is because, for given $m˙N2,const$⁠, oxygen is more easily released from the bed during reduction under higher C.

Figure 5(b) shows that δ increases monotonously in the whole range of κ while the slope increases first and then declines. The maximum η_th exists against κ corresponding to the optimal specific duration κ_opt. It implies that efficient solar-to-fuel conversion does not necessarily require large δ. The optimal efficiency η_th,max is affected by two confronting effects of the sweep gas flowrate $m˙N2,const$⁠: (1) low partial pressure of oxygen is achieved by introducing the flow of sweeping gas which can in turn promote the reduction extent of ceria and (2) gas sweeping causes heat losses which can only be partially recovered by a gas-phase heat exchanger. The optimal efficiency η_th,max is found at the stationary point where (∂η_th/∂κ) = 0. For example, for $m˙N2,const=0.005kg/s$ and C = 1500, η_th,max is found at δ = 0.028 which is less than half of the peak value of 0.07. In addition, with larger C, η_th,max and κ_opt are both increased. As shown in Fig. 5(b), η_th,max is only 2.4% obtained at κ_opt = 1890 s/kg for C = 1500, that is increased to 4.3% at larger κ_opt = 2600 s/kg for C = 3000.

The share of different energy fractions (F) varying with specific duration (κ) is given in Fig. 6. At the initial stage of reduction, sensible heat loss contributes the largest fraction, F_sen, accounting for 52% of the total energy input, followed by gas-phase heat recovery loss, F_gas, for 25%, and conductive loss, F_cond, for 19%. The radiative loss fraction, F_rad, is merely 3% at beginning. Significant change of the share is observed as κ increases to 10⁴ s/kg. F_sen shrinks to only 10% of its initial value. Meanwhile, F_gas and F_rad increase up to 50% and 27% rapidly and become the most major fractions. F_cond is less significant, accounting for 14%. It shows that the share approaches to the steady state over a certain time, i.e., κ > 2 × 10⁴ s/kg. More than 96% of the total energy input is consumed by the heat losses from gas-phase transport, cavity radiation, and conduction. F_gas, F_rad, and F_cond account for 52%, 30%, and 14%, respectively. F_chem represents the ratio of chemical energy stored in reduced ceria which determines the final thermal efficiency. During reduction, F_chem first increases, then drops down after reaching a peak value of 6.7%, and finally returns to zero.

Figure 7 shows the results of maximum thermal efficiencies (η_th,max) as functions of constant mass flows of N₂$(m˙N2,const)$⁠. The curve of black spots represents the distribution of η_th,max against κ neglecting the temperature constraint. The peak η_th,max, i.e., the optimal efficiency (η_th,opt1), is 12.6%, found at $m˙N2=3.5×10−4kg/s$⁠. However, η_th,opt1 is overestimated because reduction temperature has exceeded the upper limit for ceria material (2100 K). The second curve of spots shows results with the temperature constraint, i.e., T_max < 2100 K. The optimal efficiency (η_th,opt2) is reduced to 11.1% found at $m˙N2=8.2×10−4kg/s$⁠. The gap between two curves reflects the influence by the temperature constraint. In the actual experiments, the reactor should be operated below the temperature limit to avoid phase change of ceria. If the time it takes to reach temperature limit is shorter than the optimal operating duration, the optimal thermal efficiency will not be achieved. However, with increased $m˙N2$⁠, the two curves gradually approach and finally coincide. The overlap threshold is found at $m˙N2=2.5×10−3kg/s$⁠. For any higher sweep gas mass flowrates, the peak temperature in the reactor within the optimal operating duration will not exceed the limit of 2100 K.

The whole domain in Fig. 7 can be separated into three phases: the logarithmic phase, joint phase, and linear phase. For instance, Eq. (19) gives the piecewise fitting function according to the curve with the temperature constraint in Fig. 6. Coefficients are given in Table 2.

For the logarithmic phase, the gas-phase heat loss (Q_gas) is relatively low due to the small scale of $m˙N2$⁠. The optimal efficiency η_th,max is mainly limited by the high oxygen partial pressure. It is shown that $ηth,max∝lnm˙N2$ when $m˙N2<3×10−4kg/s$⁠, which can be explained by Eqs. (12) and (13): $m˙N2$ is inversely proportional to $PO2$⁠, which has a negative logarithmic relationship with η_th,max. For the linear phase, η_th,max is mainly limited by Q_gas which has become a dominant factor to thermal performance. From the figure, we can tell $ηth,max∝−m˙N2$ when $m˙N2>2.5×10−3kg/s$⁠. Equation (7) also shows that Q_gas increases with $m˙N2$ in linear relation. For the joint phase, the dual effect of sweeping gas co-exists and the combination is used.

Figure 8 shows the results of the final temperature (T_fin) and the radiative efficiency (η_rad)³ as functions of the mass flowrate of N₂$(m˙N2)$⁠. The curves in hollowed circles and squares show the results without the temperature constraint, T_fin reaches up to 2380 K when $m˙N2=4×10−5kg/s$ corresponding to the lowest η_rad of 67.6%. After that, T_fin decreases gradually while η_rad keeps increasing with increasing $m˙N2$⁠. The curves in solid spots and pluses represent the results with the temperature constraint. The radiative absorption efficiency η_rad is improved compared to that without the temperature constraint when $m˙N2<2.5×10−3kg/s$⁠. After that, similar to Fig. 7, the two sets of data coincide. The value of η_rad exceeds beyond 75% over the whole range of $m˙N2$⁠. The minimum value equals 77.8% found at the overlap point.

Figure 9 shows the results of the energy fractions (F) corresponding to the maximum thermal efficiency (η_th,max) in different constant mass flowrates of N₂$(m˙N2,const)$⁠. F_sen is the major fraction among all when $m˙N2<4×10−5kg/s$⁠, accounting for more than 50%, followed by F_rad = 18.6%, F_cond = 16.3%, F_chem = 13.6%, and F_gas = 12.4%. With increased $m˙N2$⁠, F_sen and F_gas vary significantly while the change of other energy fractions keeps in a relatively small range. The sensible fraction F_sen falls below 25.0% lower than F_gas as $m˙N2$ increases to 2.5 × 10⁻³ kg/s. Meanwhile, F_rad slightly increases and approaches to the peak value of 22.2%. The conduction fraction F_cond remains between 15% and 16%. The chemical energy storage fraction F_chem is reduced to 13.7% making it the smallest fraction. In the range of $m˙N2>2.5×10−3kg/s$⁠, F_gas becomes the dominant fraction which eventually exceeds 50%. Fractions F_cond, F_sen, and F_rad are similar, between 12% and 18%. By contrast, F_chem decreases rapidly and drops to 3.4% when $m˙N2=6.7×10−3kg/s$⁠.

In this section, we go further to discuss the optimization issues using variable N₂ flows $(m˙N2,var)$⁠. The mass flowrate of N₂ is under transient control during thermal reduction. GA is used for searching the globally optimal solutions.

Figures 11 and 12 show the results of temperature (T), oxygen partial pressure $(PO2)$⁠, nonstoichiometry (δ), and thermal efficiency (η_th), as the functions of the specific duration (κ) in the variable N₂ flow $(m˙N2,var)$ according to Eq. (20). As shown in Fig. 11, the reduction process is recognized into two stages: (1) T increases at a constant $PO2$ and (2) $PO2$ decreases at a constant T. For the first stage (κ < 910 s/kg), the sweeping flow is neglectable, thus air inside the cavity is almost stagnant. As generated O₂ is accumulated, $PO2$ increases rapidly and approaches 1. Reduction is driven by increasing temperature only. The nonstoichiometry of ceria (δ) increases from the initial value 3.2 × 10⁻⁵ to 0.036 during the first stage. For the second stage (κ > 910 s/kg), as T approaches to the peak value of 2100 K, $PO2$ is reduced by increasing $m˙N2$⁠. Therefore, δ is further increased under a quasi-isothermal condition, finally reaching 0.109 at T = 2100 K and $PO2=2.3×10−3$⁠.

Figure 12 shows the curves of thermal efficiency (η_th) and sweeping flowrate $(m˙N2)$ over specific durations (κ). The sweeping flowrate keeps lower than10⁻⁶ kg/s during the first stage, then rapidly increases to 2.8 × 10⁻³ kg/s at the end of the second stage. The thermal efficiency keeps increasing with κ which eventually reaches the optimal value, 12.9%. Two peaks of η_th during thermal reduction were observed. One is found in the mid-to-late first stage, ∼600 s/kg, and the other at the beginning of the second stage, ∼950 s/kg. It implies that the principle of the optimal strategy for ceria reduction in a fixed-bed model is first to increase T rapidly and then decrease $PO2$ rapidly when approaching to T_max.

Figure 13 compares two optimal cases with respect to constant and variable N₂ flows by pie diagrams. Thermal performance is improved based on the strategy of variable flows compared to that of constant flows. The gas-phase heat loss fraction (F_gas) is reduced from 7.8% to 5.5% by 29.5%. As a result, the chemical storage fraction (F_chem) increases from 17.0% to 19.3% by 13.5%.

In this paper, a transient fixed-bed model was proposed for the solar-driven ceria-based thermal reduction. Under the baseline conditions, the N₂ feeding strategies of constant and variable flows were discussed by PAM and genetic algorithm. The main conclusions are as follows:

The extreme point of the thermal efficiency corresponds to the optimal reduction duration for a given constant N₂ flowrate. Without the temperature constraint, the maximum thermal efficiency reaches 12.6% in the optimal N₂ flowrate of 3.5 × 10⁻⁴ kg/s. By contrast, the maximum efficiency drops to 11.1% with the temperature constraint, i.e., T < 2100 K. The optimal N₂ flowrate is 8.2 × 10⁻⁴ kg/s.

The flowrate was optimized with genetic algorithm for the feeding strategy of variable N₂ flowrate. Compared to the optimal case of constant N₂ flow, the thermal efficiency was increased by 16.2%. The optimal strategy involves the two steps of isobaric heating and isothermal compression.

This work was supported by the National Science Foundation of China (No. 51736006) and the Scientific Research Foundation of Graduate School of Southeast University.

m =

mass, kg

p =

absolute pressure, Pa

C =

concentration ratio

F =

energy fraction

H =

molar thermal enthalpy, J/mol

I =

direct normal irradiation, W/m²

M =

molar mass, kg/mol

P =

relative partial pressure

Q =

thermal energy, J

T =

temperature, K

$m˙$ =

mass flowrate, kg/s

x =

vector of logarithmic variables, ln(1/s)

$Q˙$ =

power, W

c_p =

specific heat at constant pressure, J/(kg K)

F_C =

conductive heat loss fraction

Δκ =

specific time-step, s/kg

Δt =

time-step, s

Δδ =

change in nonstoichiometry

$ΔhO2∘$ =

reaction molar enthalpy of ceria, J/mol(O₂)

$ΔsO2∘$ =

reaction molar entropy of ceria, J/(mol(O₂) K)

ΔH =

enthalpy difference, J

δ =

nonstoichiometry

ɛ =

heat recovery effectiveness, %

η_th =

thermal efficiency, %

κ =

specific reaction duration, s/kg

τ =

reaction duration, s

chem =

chemical

cond =

conductive

cons =

constant

fin =

finishing

gas =

gas phase

int =

initial

max =

maximum

opt =

optimal

ox =

oxidation

rad =

radiative

red =

reduction

ref =

reference

sen =

sensitive

solar =

solar energy

var =

variable

∞ =

ambient