Fuel Reactor CFD Multiscale Modelling in Syngas-Based Chemical Looping Combustion with Ilmenite

: As global power generation is currently relying on fossil fuel-based power plants, more anthropogenic CO 2 is being released into the atmosphere. During the transition period to alternative energy sources, carbon capture and storage seems to be a promising solution. Chemical-looping combustion (CLC) is an energy conversion technology designed for combustion of fossil fuel with advantageous carbon capture capabilities. In this work, a 1D computational ﬂuid dynamics (CFD) multiscale model was developed to study the reduction step in a syngas-based CLC system and was validated using literature data ( R = 0.99). In order to investigate mass transfer effects, ﬂow rate and particle dimension studies were carried out. Sharper mass transfer rates were seen at lower ﬂow rates and smaller granule sizes due to suppression of diffusion limitations. In addition, a 3D CFD particle model was developed to investigate in depth the reduction within an ilmenite particle, with focus on heat transfer effects. Minor differences of 1 K were seen when comparing temperature changes predicted by the two models during the slightly exothermic reduction reaction with syngas.


Introduction
The climate is changing at an alarmingly rapid rate as a consequence of increasing concentrations of greenhouse gases in the atmosphere. Most anthropogenic CO 2 emissions are a result of energy production based on fossil fuels. Ideally, a switch from fossil fuelbased energy production to energy production using non-fossil and other renewable sources would considerably reduce CO 2 emissions. However, during this transition period, fossil fuels will continue to be the dominant global energy source [1].
A solution that enables fossil fuel usage, while mitigating CO 2 emissions, is carbon capture and storage (CCS). CCS is the process through which CO 2 is captured at large industrial process facilities, transported to a storage site, and injected into rock formations for permanent storage. Conventional CO 2 separation techniques from diluted streams lead to significant additional energy costs and decrease in overall power efficiency [2]. This drawback can be overcome by chemical-looping combustion (CLC) [3].
CLC enables fuel combustion capable of producing a pure CO 2 exhaust stream at reduced energy penalties, achieved by avoiding direct contact between air and fuel, essentially circumventing the need to separate CO 2 and N 2 from the flue gas [4]. In CLC, a solid oxygen carrier (OC) (metal/metal oxide) acts as an intermediate between alternating oxidation and reduction reactions taking place in different reactors [5]. During oxidation, the OC reacts with air and forms an oxidized solid, producing heat, while during reduction, the oxidized OC reacts with fuel producing CO 2 and water. Following separation of water by condensation, a pure stream of CO 2 is obtained. The same amount of heat is released from CLC as from traditional combustion [6].
Conventionally, CLC takes place in an interconnected system of fluidized beds. The two reactors are an air reactor (a riser), where oxidation takes place, and a fuel reactor (a bubbling fluidized bed), where reduction takes place [7]. Drawbacks of the system are related to OC transportation due to additional energy requirements, as well as a challenging gas-particle separation at the high pressures and temperatures of the system [8].
A different approach on CLC was proposed by Noorman et al. [9], where the system uses a fixed bed reactor configuration with stationary OC. The reactor is operated dynamically by switching the gas feed streams and the OC is periodically exposed to oxidizing or reducing conditions. The main benefit to such a concept is circumventing the need for gas-particle separation due to the stationary nature of the OC. In order to mimic a continuous process, two separate parallel reactors are required [10].
This paper presents numerical modelling and simulation work supported by experimental validation with the goal of investigating the process dynamics of a fuel reactor taking part in a syngas-based CLC system. Figure 1 shows a graphical representation of the fuel reactor, as well as the three different fuel gases considered in simulations.
Energies 2021, 14, x FOR PEER REVIEW 2 of 19 Conventionally, CLC takes place in an interconnected system of fluidized beds. The two reactors are an air reactor (a riser), where oxidation takes place, and a fuel reactor (a bubbling fluidized bed), where reduction takes place [7]. Drawbacks of the system are related to OC transportation due to additional energy requirements, as well as a challenging gas-particle separation at the high pressures and temperatures of the system [8].
A different approach on CLC was proposed by Noorman et al. [9], where the system uses a fixed bed reactor configuration with stationary OC. The reactor is operated dynamically by switching the gas feed streams and the OC is periodically exposed to oxidizing or reducing conditions. The main benefit to such a concept is circumventing the need for gas-particle separation due to the stationary nature of the OC. In order to mimic a continuous process, two separate parallel reactors are required [10].
This paper presents numerical modelling and simulation work supported by experimental validation with the goal of investigating the process dynamics of a fuel reactor taking part in a syngas-based CLC system. Figure 1 shows a graphical representation of the fuel reactor, as well as the three different fuel gases considered in simulations. The novelty of the paper lies in the modelling work and results. A multiscale 1D model was developed to simulate the reduction step taking place in a CLC process with ilmenite, a titanium-iron oxide, as the oxygen carrier and was validated using breakthrough data published by Gallucci et al. [11]. Modelling work was simplified by only considering the iron component of the OC, leading to the following reduction reactions: Subsequently, a study of the system dynamics was carried out to examine the effects of varying flow rates, as well as different particle dimensions. Furthermore, a 3D CFD model was developed to study the behavior of ilmenite particles during the reduction phase of a CLC process, with focus on heat transfer effects.

Materials and Methods
In order to accurately describe the behavior of a CLC system, it is crucial to understand the phenomena occurring during the reduction step. To this end, complex CFD models accounting for mass, energy, and momentum transport were developed using COMSOL Multiphysics. The novelty of the paper lies in the modelling work and results. A multiscale 1D model was developed to simulate the reduction step taking place in a CLC process with ilmenite, a titanium-iron oxide, as the oxygen carrier and was validated using breakthrough data published by Gallucci et al. [11]. Modelling work was simplified by only considering the iron component of the OC, leading to the following reduction reactions: Subsequently, a study of the system dynamics was carried out to examine the effects of varying flow rates, as well as different particle dimensions. Furthermore, a 3D CFD model was developed to study the behavior of ilmenite particles during the reduction phase of a CLC process, with focus on heat transfer effects.

Materials and Methods
In order to accurately describe the behavior of a CLC system, it is crucial to understand the phenomena occurring during the reduction step. To this end, complex CFD models accounting for mass, energy, and momentum transport were developed using COMSOL Multiphysics. A 1D multiscale model of a fuel reactor was developed to investigate combustion with syngas. The model accounts for macroscale interactions in the packed bed, represented by the bed length, a 1D space dimension, as well as interactions at the microscale within A 1D multiscale model of a fuel reactor was developed to investigate combustion with syngas. The model accounts for macroscale interactions in the packed bed, represented by the bed length, a 1D space dimension, as well as interactions at the microscale within the particles, represented by the particle radius, and a 1D space dimension. A graphical representation of the model concept can be seen in Figure 2. Moreover, the packed bed domain was considered homogenous and was modelled as porous media. Mass, momentum, and heat transfer equations were coupled together with reduction kinetics and solved numerically in space and time using the finite element method. Three different cases were considered for the gas feed composition: H2, CO, and syngas (Table 1). The model was validated using H2 and CO breakthrough curves published by Gallucci et al. [11]. The experimental measurements were taken using a high temperature resistant steel packed bed reactor and a corresponding test rig. The lab-scale reactor could be operated at conditions of up to 1150 °C and 10 bar. Gas pre-heating to inlet temperature was carried out using ovens, while no additional heating systems were present in the reactor. By using 20 thermocouples installed at different axial positions inside the bed, as well as a mass spectrometry analysis at the reactor outlet, accurate measurement of the bed temperature and gas composition at the outlet was possible. To enable validation, the same operating conditions were considered in the simulations (Table 2). Reaction kinetics published by Parker [12] were used to simulate reduction with H2 and CO. An additional first order reaction rate was added to account for the carbon deposition observed to be occurring in the experimental measurements. Model predicted concentration profiles of the gas and solid, as well as bed temperature profiles, were investigated. Furthermore, sensitivity studies regarding flow rate and particle dimension were performed for reduction with syngas. As the multiscale model had limitations regarding particle temperature behavior and could only predict an average particle temperature, a 3D CFD model was developed to better simulate the reduction reaction taking place within an ilmenite particle. The fluid domain was represented with a cube geometry, while the ilmenite particle domain was represented using a sphere ( Figure 3). It was necessary to simulate in three dimensions of space due to the complexity of forced convection over a particle. The model accounted for Moreover, the packed bed domain was considered homogenous and was modelled as porous media. Mass, momentum, and heat transfer equations were coupled together with reduction kinetics and solved numerically in space and time using the finite element method. Three different cases were considered for the gas feed composition: H 2 , CO, and syngas (Table 1). The model was validated using H 2 and CO breakthrough curves published by Gallucci et al. [11]. The experimental measurements were taken using a high temperature resistant steel packed bed reactor and a corresponding test rig. The lab-scale reactor could be operated at conditions of up to 1150 • C and 10 bar. Gas pre-heating to inlet temperature was carried out using ovens, while no additional heating systems were present in the reactor. By using 20 thermocouples installed at different axial positions inside the bed, as well as a mass spectrometry analysis at the reactor outlet, accurate measurement of the bed temperature and gas composition at the outlet was possible. To enable validation, the same operating conditions were considered in the simulations ( Table 2). Reaction kinetics published by Parker [12] were used to simulate reduction with H 2 and CO. An additional first order reaction rate was added to account for the carbon deposition observed to be occurring in the experimental measurements. Model predicted concentration profiles of the gas and solid, as well as bed temperature profiles, were investigated. Furthermore, sensitivity studies regarding flow rate and particle dimension were performed for reduction with syngas. As the multiscale model had limitations regarding particle temperature behavior and could only predict an average particle temperature, a 3D CFD model was developed to better simulate the reduction reaction taking place within an ilmenite particle. The fluid domain was represented with a cube geometry, while the ilmenite particle domain was represented using a sphere ( Figure 3). It was necessary to simulate in three dimensions of space due to the complexity of forced convection over a particle. The model accounted for mass, momentum, and heat transport. The flow regime was established to be laminar based on the calculated Reynolds number of 320. Fluid flow and heat transfer were initially simulated in stationary state, after which the solutions from the stationary study were passed on as initial values to be used in the transient study. During reduction, grain dimensions and particle porosity were assumed constant, meaning no change with solid composition in effective diffusivity of the gas within the particle. In addition to higher complexity in temperature prediction, momentum predictions were also enhanced, by describing the laminar fluid flow surrounding the particle and its effects over mass transport. The computational domain for fluid flow was considered as the void between the particles. The dimension of the fluid flow region was calculated by assuming evenly distributed pellets in the reactor bed. While the 3D single particle model had certain constraints, such as not being able to capture the underlying complexities of an entire packed bed of particles, the scope of the 3D single particle model was the in-depth analysis of the intraparticle phenomena occurring during the reduction reaction, with strong emphasis on heat transport related effects. mass, momentum, and heat transport. The flow regime was established to be laminar based on the calculated Reynolds number of 320. Fluid flow and heat transfer were initially simulated in stationary state, after which the solutions from the stationary study were passed on as initial values to be used in the transient study. During reduction, grain dimensions and particle porosity were assumed constant, meaning no change with solid composition in effective diffusivity of the gas within the particle. In addition to higher complexity in temperature prediction, momentum predictions were also enhanced, by describing the laminar fluid flow surrounding the particle and its effects over mass transport. The computational domain for fluid flow was considered as the void between the particles. The dimension of the fluid flow region was calculated by assuming evenly distributed pellets in the reactor bed. While the 3D single particle model had certain constraints, such as not being able to capture the underlying complexities of an entire packed bed of particles, the scope of the 3D single particle model was the in-depth analysis of the intraparticle phenomena occurring during the reduction reaction, with strong emphasis on heat transport related effects.  For both 1D multiscale and 3D particle models, the mesh was sufficiently refined until mesh convergence was achieved. Figure 3 displays the high quality of the mesh in terms of element skewness and the increased density of elements discretizing the particle domain. Additionally, the tolerance factor was tightened accordingly to guarantee solution accuracy. For both 1D multiscale and 3D particle models, the mesh was sufficiently refined until mesh convergence was achieved. Figure 3 displays the high quality of the mesh in terms of element skewness and the increased density of elements discretizing the particle domain. Additionally, the tolerance factor was tightened accordingly to guarantee solution accuracy.

Model Assumptions
Model assumptions are presented below: • the packed bed domain was considered homogenous; • the solid iron-based OC particles were considered stationary; • the reduction reaction occurred within the OC particles; • the OC particles were identical, with a spherical shape and constant radius;  (3)) in both models was: where ρ was fluid density and u was velocity.

Ergun's Equation
The packed bed in the fuel reactor model was assumed to be a uniform porous domain filled with ilmenite OC granules. The distribution of pressure along the reactor length was solved with the differential form of Ergun's equation (Equation (4)) [13]: where ε b was bed voidage, D pe was pellet diameter and µ was fluid viscosity.

Navier-Stokes and Brinkman Equations
The particle model used the Navier-Stokes equations (Equation (5)) together with the continuity equation (Equation (3)) to solve laminar fluid flow surrounding the granule: where τ was the viscous stress tensor and was calculated with Equation (6): The Brinkman equations (Equation (7)) were used to describe porous media flow in the laminar regime within the ilmenite particle: where ε p was particle porosity and κ was permeability.

Species Transport
The fuel reactor model solved a convection-diffusion equation (Equation (8)) to represent mass transport in the bed: where c was the concentration in the bulk and D L was the axial dispersion coefficient estimated with a correlation for packed beds (Equation (9)), published by Wakao and Funazkri [14], based on the Reynolds, Peclet, and Schimdt numbers: Mass transfer within the particle in the multiscale model was represented by diffusion of gas components as a result of a concentration gradient, as well as the reduction reaction between the OC solid and fuel. A diffusion-reduction reaction equation (Equation (10)) was solved along a dimensionless radial coordinate ranging from the center to the particle surface for each of the gas components: and for the solid OC components (Equation (11)): where r was a dimensionless radial coordinate, r p was the particle radius, N was the number of granules per unit volume of bed, R was the reduction reaction rate described with kinetics published by Parker (Equations (21) and (22)) [12], and D p e f f was the effective diffusion coefficient within the particle calculated using the Millington and Quirk model (Equation (12)): The coupling between macroscale and microscale dimensions was achieved by assuming film resistance at the fluid-pellet interface, calculating the surface flux with Equation (13): where h D,i was the mass transfer coefficient calculated with Equation (14): For the 3D particle model, mass transport in the fluid domain was described with a convection-diffusion equation (Equation (15)): Within the particle, a diffusion-reduction reaction was used to describe mass transfer of the gas components (Equation (16)): and for the solid OC components (Equation (17)): Intraparticle diffusion (Equation (12)) and reaction rate (Equations (21) and (22)) were calculated with the same expressions as in the multiscale model. Film resistance at the fluid-particle interface was also accounted for and calculated with Equation (14).

Energy Transport
For both multiscale and particle models, a homogenized approach was chosen to describe heat transport by assuming the local thermal equilibrium hypothesis [15], meaning local fluid temperature was considered identical to porous media temperature and only one temperature variable needed to be computed at each spatial location in the domain. A single heat transfer equation was solved (Equation (18)): The effective material properties accounting for both fluid and porous domain properties were defined through volumetric averaging. As such, for effective heat capacity (Equation (19)): ρC and for effective thermal conductivity (Equation (20)): The external term Q R was either a heat source or sink term depending on the enthalpy of reaction and the gas composition chosen for reduction with the ilmenite OC.

Reduction Reaction Kinetics
For both multiscale and particle models, the kinetic model used to calculate the reduction reactions with the ilmenite OC was published by Abad et al. [16]. The model assumed process velocity was subject to reaction rate, neglecting mass transfer effects related to pellet-fluid film resistance and intraparticle diffusion, both of which were accurately described in the CFD models. The following equations, subsequently adjusted by Parker et al. [12], were used in this work for defining reduction with H 2 (Equation (21)): and reduction with CO (Equation (20)): Moreover, to highlight the carbon deposition effect observed in the experimental measurements by Gallucci et al. [11], a first order reaction rate as a function of CO concentration was added to the models (Equation (23)): where k CO,d was the first order rate constant determined using published experimental data [11].

Model Parameters and Operating Conditions
In order to enable model validation, the parameters used in the models were taken from a paper in which Gallucci et al. [11] published experimental measurements of a syngasbased CLC process with ilmenite as OC. Process, operating, and geometry parameters are shown in Table 2.

Results and Discussion
The model simulating iron-based OC reduction with H 2 was validated using H 2 breakthrough curves published by Gallucci et al. [11]. Figure 4 shows the comparison between experimental and model predicted H 2 breakthrough curves.

Model Validation-Reduction with H2
The model simulating iron-based OC reduction with H2 was validated using H2 breakthrough curves published by Gallucci et al. [11]. Figure 4 shows the comparison between experimental and model predicted H2 breakthrough curves. It can be noticed that the experimental measurements showed an increase in H2 concentration over the feed concentration when steady state was reached. Gallucci et al. [11] reasoned that a steam-iron process was occurring during the experiments and led to some hydrogen being produced during reduction. While the steam-iron process was not accounted for in the model, a good agreement was seen between experimental and model predicted breakthrough curves when calculating the Pearson coefficient correlation ( ) with a value of 0.9995. Figure 5 shows H2 bulk concentrations and concentrations in the center of the ironbased OC particles represented as a function of reactor axial position at various simulation times.
The difference between H2 concentration in the bulk and center of the particle was attributed to mass transfer effects, such as film resistance at the fluid-particle interface and intraparticle diffusion. Furthermore, the iron-based OC was reduced as the H2 gas diffused further within the particle. Figure 6 presents the change in concentration of Fe2O3 and FeO within four ironbased OC particles selected at different axial positions of 0.2 m, 0.4 m, 0.6 m, 0.8 m, and at a simulation time of 700 s. The normalized concentrations of Fe2O3 and FeO were represented as function of normalized particle radius. The normalized particle radius was defined as the ratio of the spatial radial coordinate in the particle to the particle radius, meaning that the particle modelling domain ranged from 0 (centre) to 1 (particle perimeter). The normalized concentrations of the oxides were calculated by dividing Fe2O3 concentration at simulation time by the initial concentration of Fe2O3 and by dividing FeO concentration at simulation time by the final concentration of FeO. It can be noticed that the experimental measurements showed an increase in H 2 concentration over the feed concentration when steady state was reached. Gallucci et al. [11] reasoned that a steam-iron process was occurring during the experiments and led to some hydrogen being produced during reduction. While the steam-iron process was not accounted for in the model, a good agreement was seen between experimental and model predicted breakthrough curves when calculating the Pearson coefficient correlation (R) with a value of 0.9995. Figure 5 shows H 2 bulk concentrations and concentrations in the center of the iron-based OC particles represented as a function of reactor axial position at various simulation times.  The difference between H 2 concentration in the bulk and center of the particle was attributed to mass transfer effects, such as film resistance at the fluid-particle interface and intraparticle diffusion. Furthermore, the iron-based OC was reduced as the H 2 gas diffused further within the particle.  3 and FeO were represented as function of normalized particle radius. The normalized particle radius was defined as the ratio of the spatial radial coordinate in the particle to the particle radius, meaning that the particle modelling domain ranged from 0 (centre) to 1 (particle perimeter). The normalized concentrations of the oxides were calculated by dividing   It can be noticed that Fe2O3 was fully converted to FeO and reduction had ended within the iron OC particle at 0.2 m, while the reduction reaction was still ongoing for the OC particles situated at 0.4 m and 0.6 m. In the case of the final OC particle at 0.8 m, the H2 gas had not reached the particle thus far and reduction had not occurred. The changes seen in Figure 6 can be correlated with the H2 concentration presented in Figure 5 at simulation time 700 s (green line).
Temperature profiles, predicted from the multiscale model during the reduction reaction of ilmenite with H2, are shown in Figure 7. It can be noticed that Fe 2 O 3 was fully converted to FeO and reduction had ended within the iron OC particle at 0.2 m, while the reduction reaction was still ongoing for the OC particles situated at 0.4 m and 0.6 m. In the case of the final OC particle at 0.8 m, the H 2 gas had not reached the particle thus far and reduction had not occurred. The changes seen in Figure 6 can be correlated with the H 2 concentration presented in Figure 5 at simulation time 700 s (green line).
Temperature profiles, predicted from the multiscale model during the reduction reaction of ilmenite with H 2 , are shown in Figure 7.
In their experimental measurements, Gallucci et al. [11] used G11 ilmenite particles with a weight percentage of 35% iron oxide and 65% titanium (see Table 2) and reported a rise in temperature during the reduction experiments with H 2 , meaning an exothermic reaction took place. This work considered only the iron component of the ilmenite OC to describe the reduction reactions (see Equations (1) and (2)) and the redox pair of FeO and Fe 2 O 3 would have been normally used to describe an endothermic reaction with H 2 ; however, in order to match the experimental measurements, the enthalpies of reaction of the redox pair of FeOTiO 4 and Fe 2 TiO 5 were selected to calculate change in temperature. Therefore, it can be noticed that the model predicted a temperature increase along the reactor length, from an initial temperature of 873 K, as a result of the exothermic reduction reaction of ilmenite with H 2 . In their experimental measurements, Gallucci et al. [11] used G11 ilmenite particles with a weight percentage of 35% iron oxide and 65% titanium (see Table 2) and reported a rise in temperature during the reduction experiments with H2, meaning an exothermic reaction took place. This work considered only the iron component of the ilmenite OC to describe the reduction reactions (see Equations (1) and (2)) and the redox pair of FeO and Fe2O3 would have been normally used to describe an endothermic reaction with H2; however, in order to match the experimental measurements, the enthalpies of reaction of the redox pair of FeOTiO4 and Fe2TiO5 were selected to calculate change in temperature. Therefore, it can be noticed that the model predicted a temperature increase along the reactor length, from an initial temperature of 873 K, as a result of the exothermic reduction reaction of ilmenite with H2.

Model Validation-Reduction with CO
The model describing iron-based OC reduction with CO was validated using CO breakthrough curves published by Gallucci et al. [11]. Figure 8 shows the comparison between experimental and model predicted CO breakthrough curves.

Model Validation-Reduction with CO
The model describing iron-based OC reduction with CO was validated using CO breakthrough curves published by Gallucci et al. [11]. Figure 8 shows the comparison between experimental and model predicted CO breakthrough curves. In their experimental measurements, Gallucci et al. [11] used G11 ilmenite particles with a weight percentage of 35% iron oxide and 65% titanium (see Table 2) and reported a rise in temperature during the reduction experiments with H2, meaning an exothermic reaction took place. This work considered only the iron component of the ilmenite OC to describe the reduction reactions (see Equations (1) and (2)) and the redox pair of FeO and Fe2O3 would have been normally used to describe an endothermic reaction with H2; however, in order to match the experimental measurements, the enthalpies of reaction of the redox pair of FeOTiO4 and Fe2TiO5 were selected to calculate change in temperature. Therefore, it can be noticed that the model predicted a temperature increase along the reactor length, from an initial temperature of 873 K, as a result of the exothermic reduction reaction of ilmenite with H2.

Model Validation-Reduction with CO
The model describing iron-based OC reduction with CO was validated using CO breakthrough curves published by Gallucci et al. [11]. Figure 8 shows the comparison between experimental and model predicted CO breakthrough curves. The feed concentration of CO was 30 mol% (see Table 1, case 2), while the concentration at the outlet was 23 mol%. The supplementary consumption of CO was measured experimentally and was attributed to carbon deposition during the reduction reaction. Therefore, an additional first order reaction rate was added to the model to account for carbon deposition (Equation (23)). The slight misalignment of the lower section can be attributed to an overestimation of axial dispersion by the model. The reason for the deviation seen in the upper sections of the breakthrough curves was assumed to be due to the mass transfer limitations, unaccounted for by the first order reaction used to describe carbon deposition. Normally, operation of the CLC reactor should be performed in such a manner to suppress carbon deposition; for instance, adding steam to the syngas composition. The calculated Pearson correlation coefficient (R) of 0.9971 showed a good agreement between experimental and model predicted breakthrough curves. Figure 9 presents CO concentrations in the bulk and in the center of the iron-based OC particles represented as a function of reactor axial position at various simulation times. experimentally and was attributed to carbon deposition during the reduction reaction. Therefore, an additional first order reaction rate was added to the model to account for carbon deposition (Equation (23)). The slight misalignment of the lower section can be attributed to an overestimation of axial dispersion by the model. The reason for the deviation seen in the upper sections of the breakthrough curves was assumed to be due to the mass transfer limitations, unaccounted for by the first order reaction used to describe carbon deposition. Normally, operation of the CLC reactor should be performed in such a manner to suppress carbon deposition; for instance, adding steam to the syngas composition. The calculated Pearson correlation coefficient ( ) of 0.9971 showed a good agreement between experimental and model predicted breakthrough curves. Figure 9 presents CO concentrations in the bulk and in the center of the iron-based OC particles represented as a function of reactor axial position at various simulation times. The concentrations in the bulk and at the center of the particle were close to overlapping, with a correlation coefficient ( ) of 1, meaning the reaction rate was slow and the gas diffused quicker within the particle than reduction took place. Moreover, by comparing with the concentration profiles in Figure 5, H2 showed a higher reaction rate than CO.
Initial reactor temperature during reduction with CO was 873 K. A slight decline in temperature can be seen in Figure 10 along the reactor length at various simulation times, as a result of the endothermic reduction reaction of ilmenite with CO. The concentrations in the bulk and at the center of the particle were close to overlapping, with a correlation coefficient (R) of 1, meaning the reaction rate was slow and the gas diffused quicker within the particle than reduction took place. Moreover, by comparing with the concentration profiles in Figure 5, H 2 showed a higher reaction rate than CO.
Initial reactor temperature during reduction with CO was 873 K. A slight decline in temperature can be seen in Figure 10 along the reactor length at various simulation times, as a result of the endothermic reduction reaction of ilmenite with CO.

Reduction with Syngas
As a result of validation of the models describing iron-based OC reduction with H2 and CO, a model simulating ilmenite reduction with syngas was developed. The composition of syngas considered can be seen in Table 1, case 3. Moreover, additional simulations considering flow rates with a standard deviation of 20 NL·min −1 relative to the initial flow rate (see Table 2), as well as particle dimensions with a standard deviation of 2.76

Reduction with Syngas
As a result of validation of the models describing iron-based OC reduction with H 2 and CO, a model simulating ilmenite reduction with syngas was developed. The composition of syngas considered can be seen in Table 1, case 3. Moreover, additional simulations considering flow rates with a standard deviation of 20 NL·min −1 relative to the initial flow rate (see Table 2), as well as particle dimensions with a standard deviation of 2.76 mm from the initial particle dimension (Table 2), were carried out.
It can be noticed that breakthrough of H 2 and CO, the moment at which the values of H 2 and CO concentrations in the effluent remained unchanged due to the oxygen from the OC being entirely consumed and the reduction reactions ending (see Figure 11), occurred in much less time after reduction with syngas when compared to the breakthrough times following reduction with the gases separately (Figures 4 and 8). The same flow rate and mass of OC was considered in all three cases; however, the presence of both H 2 and CO consuming the oxygen provided by the OC, as well as higher concentrations, naturally led to a quicker breakthrough of the gases. When looking at the different flow rates, breakthrough was achieved the quickest at the highest flow rate due to H 2 and CO traveling at a higher velocity and consuming the available oxygen in less time. On the other hand, the largest breakthrough time was seen in the case of the lowest flow rate.   A similar effect can also be observed when investigating the different CO concentrations due to carbon deposition at steady state, after breakthrough. A short residence time translates into less carbon being deposited during reduction with CO.
In order to investigate mass transfer efficiency, normalized concentration was plotted against normalized time in Figure 12. Normalized concentration for H 2 was calculated by dividing concentration at simulation time with inlet concentration, while for CO it was calculated by dividing concentration with concentration at steady state, after breakthrough (due to carbon deposition). The normalized time was calculated by dividing the simulation time with the time corresponding to a breakthrough point of 0.5.  Assessment of the shape of the breakthrough curves showed steeper profiles at all flow rates for H2, rather than for CO, meaning faster kinetics, a smaller mass transfer zone, and a higher reaction rate when reducing iron-based OC with H2, corresponding with findings published by Ortiz et al. [17] on redox kinetics regarding ilmenite granules used in CLC. Finally, with increasing flow rate, the breakthrough profiles lost sharpness, which was expected behavior.
The effects of particle dimension over mass transfer rate can be seen in Figure 13. The mass of ilmenite OC was identical for all cases; as such, the breakthrough times were also identical. In order to investigate mass transfer efficiency, normalized concentration at the outlet was plotted against normalized time. The sharpest profile for both CO and H2 was seen when considering the smallest particle dimension due to the shorter diffusion path length inside the particle, leading to suppressing diffusion limitations and maximizing the reduction reaction inside the OC. This then indicates that mass transfer efficiency increased with a decrease in particle size. Assessment of the shape of the breakthrough curves showed steeper profiles at all flow rates for H 2 , rather than for CO, meaning faster kinetics, a smaller mass transfer zone, and a higher reaction rate when reducing iron-based OC with H 2 , corresponding with findings published by Ortiz et al. [17] on redox kinetics regarding ilmenite granules used in CLC. Finally, with increasing flow rate, the breakthrough profiles lost sharpness, which was expected behavior.
The effects of particle dimension over mass transfer rate can be seen in Figure 13.   The mass of ilmenite OC was identical for all cases; as such, the breakthrough times were also identical. In order to investigate mass transfer efficiency, normalized concentration at the outlet was plotted against normalized time. The sharpest profile for both CO and H 2 was seen when considering the smallest particle dimension due to the shorter diffusion path length inside the particle, leading to suppressing diffusion limitations and maximizing the reduction reaction inside the OC. This then indicates that mass transfer efficiency increased with a decrease in particle size.    It can be noticed that the most Fe2O3 converted to FeO was within the iron OC particle at 0.2 m, while the reduction reaction was still ongoing for all the other OC particles situated at 0.4 m, 0.6 m, and 0.8 m. The concentrations of the oxides were much closer between the selected particles than when compared to reduction with H2 ( Figure 6), due to the fact that oxygen consumption during reduction with syngas was much higher. Figure 15 presents the change in temperature along the reactor length at various simulation times during reduction with syngas. There was an overall rise in temperature as the exothermic reduction with H2 released more energy than was required by the endothermic reduction with CO. It can be noticed that the most Fe 2 O 3 converted to FeO was within the iron OC particle at 0.2 m, while the reduction reaction was still ongoing for all the other OC particles situated at 0.4 m, 0.6 m, and 0.8 m. The concentrations of the oxides were much closer between the selected particles than when compared to reduction with H 2 (Figure 6), due to the fact that oxygen consumption during reduction with syngas was much higher. Figure 15 presents the change in temperature along the reactor length at various simulation times during reduction with syngas. There was an overall rise in temperature as the exothermic reduction with H 2 released more energy than was required by the endothermic reduction with CO.

Single Particle Model (3D)-Reduction with Syngas
In order to further study the behavior of ilmenite particles during the reduction phase of the CLC process, a 3D CFD transient model was developed. Fluid flow was initially solved in stationary state. Cross-sectional surfaces highlighting the velocity profiles can be seen in Figure 16.

Single Particle Model (3D)-Reduction with Syngas
In order to further study the behavior of ilmenite particles during the reduction phase of the CLC process, a 3D CFD transient model was developed. Fluid flow was initially solved in stationary state. Cross-sectional surfaces highlighting the velocity profiles can be seen in Figure 16.

Single Particle Model (3D)-Reduction with Syngas
In order to further study the behavior of ilmenite particles during the reduction phase of the CLC process, a 3D CFD transient model was developed. Fluid flow was initially solved in stationary state. Cross-sectional surfaces highlighting the velocity profiles can be seen in Figure 16. A boundary condition assuming fully developed laminar flow was considered at the inlet of the fluid domain. The slices in Figure 16 accurately described the laminar nature of the flow, with a lower velocity near the edges and inside of the particle, as well as near the edges of the fluid domain, while a higher velocity was seen in the center of the fluid sections surrounding the particle.
The evolution of the OC particle during reduction can be observed in Figure 17, as well as the interdependency between the phenomena accounted for in the model. A boundary condition assuming fully developed laminar flow was considered at the inlet of the fluid domain. The slices in Figure 16 accurately described the laminar nature of the flow, with a lower velocity near the edges and inside of the particle, as well as near the edges of the fluid domain, while a higher velocity was seen in the center of the fluid sections surrounding the particle.
The evolution of the OC particle during reduction can be observed in Figure 17, as well as the interdependency between the phenomena accounted for in the model.
The relative concentration of Fe 2 O 3 was represented by slices within the particle domain, seen in the center of the fluid domain, while the relative concentration of FeO can be observed outside of the fluid domain, also represented by slices of the particle. Moreover, the path of flow and velocity field can be visualized using streamlines with color and thickness proportional to the velocity magnitude. Finally, Figure 17 highlights temperature distribution at the outlet boundary with a temperature contour plot. Initially at simulation time of 0 s (top), only nitrogen was present in the fluid domain and Fe 2 O 3 was unconsumed. After syngas entered the fluid domain and diffused within the particle, the reduction reaction began and oxygen was consumed from the OC. At simulation time of 30 s (middle), we notice that Fe 2 O 3 concentration decreased, while FeO concentration increased proportionately, according to reaction kinetics. The effects of forced convection were also noticeable, as the heat released within the particle due to the slightly exothermic reduction reaction was transported and the temperature at the outlet boundary increased to 876 K. At 150 s (bottom), the reduction reaction was near its end and the oxygen from the OC was almost entirely consumed, as evidenced by the relative concentration of the FeO of 0.94. Furthermore, the colder inlet temperature and the slower reaction rate led to a decrease in the temperature at the outlet boundary.  The relative concentration of Fe2O3 was represented by slices within the particle domain, seen in the center of the fluid domain, while the relative concentration of FeO can be observed outside of the fluid domain, also represented by slices of the particle. Moreover, the path of flow and velocity field can be visualized using streamlines with color and thickness proportional to the velocity magnitude. Finally, Figure 17 highlights temperature distribution at the outlet boundary with a temperature contour plot. Initially at simulation time of 0 s (top), only nitrogen was present in the fluid domain and Fe2O3 was Heat transfer and fluid flow were solved together in steady state, after which the stationary solution was used as an initial value for the transient solver to avoid inconsistent initial values. Figure 18 shows an isothermal surface plot of temperature within the fluid and particle domains, as well as a contour plot of temperature at the bottom boundary of the fluid domain. a decrease in the temperature at the outlet boundary.
Heat transfer and fluid flow were solved together in steady state, after which the stationary solution was used as an initial value for the transient solver to avoid inconsistent initial values. Figure 18 shows an isothermal surface plot of temperature within the fluid and particle domains, as well as a contour plot of temperature at the bottom boundary of the fluid domain. The highest temperature change was seen at 30 s, with a maximum temperature of 878 K within the particle as a result of the exothermic reduction reaction. The released heat diffused outside of the particle and was carried away by the flow, as evidenced by the isothermal surface plot surrounding the particle. Moreover, the temperature of the bottom boundary of the fluid domain had also increased. Comparing the heat transfer results of the particle model with the temperature increase predicted by the multiscale model at 30 s of 878 K (marked with black triangle symbols in Figure 15), the difference between model predictions was only 1 K.

Conclusions
A 1D multiscale CFD packed bed reactor model within a syngas-based CLC system with ilmenite OC granules was developed using COMSOL Multiphysics. The model described the reduction step, with equations accounting for mass, heat, and momentum transfer. The model predicted breakthrough curves for reduction with H2 ( = 0.999) and CO ( = 0.997) showed a good agreement with published experimental data. Subsequently, the model was used to study flow rate and particle dimension effects over mass The highest temperature change was seen at 30 s, with a maximum temperature of 878 K within the particle as a result of the exothermic reduction reaction. The released heat diffused outside of the particle and was carried away by the flow, as evidenced by the isothermal surface plot surrounding the particle. Moreover, the temperature of the bottom boundary of the fluid domain had also increased. Comparing the heat transfer results of the particle model with the temperature increase predicted by the multiscale model at 30 s of 878 K (marked with black triangle symbols in Figure 15), the difference between model predictions was only 1 K.

Conclusions
A 1D multiscale CFD packed bed reactor model within a syngas-based CLC system with ilmenite OC granules was developed using COMSOL Multiphysics. The model described the reduction step, with equations accounting for mass, heat, and momentum transfer. The model predicted breakthrough curves for reduction with H 2 (R = 0.999) and CO (R = 0.997) showed a good agreement with published experimental data. Subsequently, the model was used to study flow rate and particle dimension effects over mass transfer rate when considering syngas as feed gas. A higher flow rate achieved breakthrough sooner, but the mass transfer rate was less sharp when compared to lower flow rates. The particle size sensitivity analysis indicated a more efficient mass transfer with decreasing particle size due to suppressing diffusion limitations and maximizing reduction inside the OC particle.
Additionally, a 3D CFD particle model was developed to enable a detailed investigation into the reduction step at the microscale level. The model described the reduction reaction using syngas within an ilmenite particle, with equations accounting for mass, heat, and momentum transfer. The individual phenomena involved, as well as the interdependency between them, were studied and discussed. As opposed to the multiscale model which considered a dynamic reactor temperature, but a constant particle temperature, the particle model's complexity led to an accurate description of heat transfer on a microscale level. However, comparison of temperature changes predicted by the models indicated a small difference of just 1 K during the exothermic reduction reaction under identical operating conditions. Future work will focus on model development in order to simulate the oxidation step, thus enabling the simulation of a full cycle of the syngas-based CLC system with ilmenite OC, followed by process optimization studies to extract useful data for process scale-up.