Numerical Analyses of Heterogeneous CLC Reaction and Transport Processes in Large Oxygen Carrier Particles

: Heterogeneous chemical looping combustion (CLC) reactions and conjugate transports in large oxygen carrier particles were numerically investigated with computational ﬂuid dynamics (CFD) approaches, in which a simpliﬁed noncatalytic reaction model was implemented for reducing intraparticle modelling computation. Volumic gas-solid reactions were treated as surface reactions based on the equivalent internal surface in the particle model. In large porous particles such as ﬁxed bed CLC reactors, the heterogeneous reactions are often limited by intraparticle diffusion. Comprehensive analyses were conducted on transports across the particle surface and their inﬂuences on reactions inside the single particles. A threshold Reynolds number of external convections was found for the enhancement of intraparticle reactions. The heterogeneous reactions, intraparticle diffusions and interstitial transports in a ﬁxed bed CLC reactor randomly packed with 597 spheres were thoroughly analysed with the same numerical approaches. Comprehensive insights of the temporal evolution and spatial distribution of scalars in the packed bed reactor were presented.


Introduction
As a novel clean combustion technology, chemical looping combustion (CLC) is attractive for the feature of inherent CO 2 capture, together with zero NOx emission and high combustion efficiency. In a CLC process, direct contact between air and fuel is avoided by dividing the conventional combustion process into two steps, performed in an air reactor and a fuel reactor, in which the oxygen carrier is oxidized and reduced, respectively. Meanwhile, the reaction absorbs or releases heat. The oxygen carrier is cycled and reused in the alternate and repetitive reduction-oxidization processes. The products of the reduction step (mainly water vapor and CO 2 ) can be condensed to remove water vapor so as to separate carbon dioxide conveniently without additional energy consumption.
In the past two decades, plenty of investigations on CLC, implemented with either fluidized bed or fixed bed reactors at laboratory and pilot scales, have been conducted [1][2][3][4][5][6]. Fluidized bed-based CLC, typically with a riser as an air reactor and a bubbling bed as a fuel reactor, is characterized by advantages such as higher heat transfer efficiency, continuous operation and stable production of hot gas. Fixed bed implementation, in which stationary oxygen carrier particles are alternately exposed to oxidation and reduction streams, was proposed later to avoid some shortcomings of its fluidized bed counterpart, such as the inevitable attrition of particles and the challenges of removing fine particles from hot stream in the conventional reactors [7].
As a novel technology concept, chemical looping combustion has been investigated experimentally and numerically to thoroughly understand the relevant reaction and transport processes which occur in the reactors of either the fluidized bed or fixed bed. Comparably, experimental studies on fixed bed CLC might be faced with more challenges due to some measurement difficulties, such as temperature and concentration distribution (as demonstrated by the axial temperature measurement by Noorman et al. [8] and the axial and radial temperature measurement by Guo et al. [6]). In general, numerical simulation can be an effective approach to overcome the difficulties faced in experimental investigations. The one-dimensional reactor model-based effective parameter method and homogeneous assumptions are usually applied to predict the temporal variation and axially spatial distribution of quantities such as temperature and species concentration during the transient fixed bed-based CLC process [9][10][11]. The Aspen Plus platform, which takes the thermodynamic equilibrium model as the kernel, has been widely used in the process simulation of thermodynamic evaluation and sensitivity analysis in chemical looping combustion and reforming systems in recent years [12,13]. On the other hand, simulations based on the three-dimensional computational fluid dynamics (CFD) method with discrete particles provide much more detail of the local flow fields and transport phenomena in fixed beds by solving the governing partial differential equations of flows and transport processes. In the past decades, a considerable number of scholars have carried out discrete particle CFD simulations, extending from detailed flow [14][15][16] to coupling heat transfer [17,18] as well as reactions [19][20][21][22].
The CLC in the fixed bed is a complicated transient process that includes heterogeneous noncatalytic reaction and species diffusion occurring inside oxygen carrier (OC) particles and the interstitial flows together with heat and mass transfer. Therefore, modelling of heterogeneous reactions in OC particles must be included when CFD-based simulation is conducted. Khakpoor et al. used a grain model (GM) based on a series of reasonable grain assumptions to determine the kinetic parameters of ilmenite reduction reaction with high accuracy [23]. In our previous work [24], the changing grain size model (CGSM) was integrated into an overset based CFD model of CLC in a random packed bed. The CGSM [25,26] assumes that the active materials of fine and nonporous grains are uniformly distributed in the porous particles, in which the grains grow or shrink in reactions following a shrinking core model (SCM). Meanwhile, the heterogeneous reaction rate and the overall particle porosity are constantly updated with the change of grain size, and this complex microscopic kinetics model makes the computational costs tremendously large and limits its application in engineering. Therefore, it is of great significance to adopt a simplified heterogeneous microscopic reaction model and improve the computational efficiency of simulation.
Actually, random packed beds with large particles (leading to a small D/d ratio) can be reasonable configurations for many applications considering their small pressure drop and enhanced radial heat removal [27], which can also be suitable for situations with exothermic reactions, such as CLC. In relatively large porous particles, species diffusion is often a limiting step for the heterogeneous reactions, as has been discussed in the works of Noorman et al. [28] and Guo & Zhu [24]. These authors have revealed that internal diffusion in large OC particles significantly affects the overall heterogeneous reaction rate, and good choices of particle properties (such as internal pore size and porosity) and particle structure (e.g., eggshells, Guo & Zhu [24]) could alleviate the limitation of internal diffusion and improve the conversion of the OC. Therefore, comprehensive and accurate numerical analyses of the effects of intraparticle diffusion characteristics in large particles on the heterogeneous reactions could effectively help to improve the gas-solid reaction performance and the utilization of the OC.
In this work, a simplified heterogeneous noncatalytic gas-solid reaction model was implemented in a commercial flow solver, attempting to reduce the computational costs related to intraparticle modelling. Numerical analyses of a large single OC particle were carried out to understand the effects of external convection and particle properties on the intraparticle diffusion characteristics and conjugated transports in the process of CLCbased CuO. The heterogeneous reactions, along with intraparticle diffusion and interstitial transports involved in CLC process in a random packed bed of 597 porous spheres, were also thoroughly analysed with the simplified intraparticle modelling approach. The com-prehensive analyses of the spatial distribution and temporal evolution of scalars helps us to understand the heterogeneous CLC reaction and transport processes in the packed bed.

The Numerical Methods and the Physical Models
In this paper, the commercial solver ANSYS Fluent, implemented with a simplified gas-solid noncatalytic reaction model, was utilized to solve the transient heterogeneous reactions in large OC particles that are involved in flow fields with conjugate transports, and all simulations were performed in parallel on a computing cluster.

Numerical Methods
In the present work, a three-dimensional heterogeneous reaction and transport processes in CLC reactor were solved numerically. The conservation equations of mass, momentum and energy for the fluid domain can be found in the Fluent theory manual [29]. To accurately solve the turbulent flow fields around particles, the k-ω SST turbulence model was adopted due to its accuracy performance for low Reynolds number flows. Moreover, a fine mesh (y+ ≤ 1) was warranted near the walls for accurate resolution.
In particle domains, we assumed an OC particle as a porous medium, and only species diffusion and heat conduction were considered. Therefore, the conservation equations of energy and species transports were used as governing equations that describe the heat and mass transfer and heterogeneous chemical reaction processes inside the particles, as shown in Equations (1) and (2).
where λ eff is the effective thermal conductivity inside porous particles, as shown in Equation (3), h is the mixture specific enthalpy, as shown in Equation (4); and S h and S m,i are the source terms of heat and species generated by heterogeneous chemical reactions, respectively, as shown in Equations (5) and (6).
The diffusive flux J i , shown in Equation (7), and Fick's law are considered in both fluid and particle domains. Inside the particle domain, only the mass diffusion equation is solved, since convection is often negligible in porous particles. In addition, some relevant works (e.g., Noorman et al. [30]) have disucssed that Knudsen diffusion in porous particles usually has more significant influence than molecular diffusion. According to the correlation by Zeng et al. [31], the mean free path of gas molecules (λ > 10d pore ) in the present work was calculated, and the diffusion resistance of gas molecules inside particles was mainly determined by the collision between molecules and the inner pore walls. Thus, the molecular diffusion resistance was ignored in the porous particle domain. As shown in Equation (8), the diffusivity terms Г i are needed for both the fluid and particle domains, which are controlled by user-defined functions (UDFs).
where D m,i is the molecular diffusivity using the correlation proposed by Fuller et al. [32] The turbulent diffusivity was calculated as µ t /Sc t,i , in which µ t and Sc t,i are turbulent viscosity and Schmidt number (0.7 for current work). D Knud,i is the Knudsen diffusion coefficient of species i, as shown in Equation (9).
where d pore is the mean pore diameter and M i is the molecular weight of species i. In this work, the transient heterogeneous chemical reaction was numerically simulated using the surface reaction of porous media in the species transport model, in which volumic heterogeneous reactions were integrated as a surface reaction on the equivalent internal surface of the porous media.
Neglecting the reverse reaction, the net heterogeneous chemical reaction rate was calculated as shown in Equation (10). The kinetic constant k s of chemical reaction was calculated with the Arrhenius equation, which is empirically related to temperature, as shown in Equation (11).
where k 0 is the pre-exponential factor, β is the temperature exponent and E act is reaction activation energy. S v is the surface-to-volume ratio of oxygen carrier, which is the ratio of the active surface area to the finite-volume cell volume. Here, the simplified heterogeneous reaction model [33] of η was implanted as shown in Equation (12), which describes the local abundance of active materials per unit volume of porous particles, and X s is the solid reactants conversion. The one-step oxidation and reduction reaction of copper-based CLC is shown in Equations (13) and (14), and the detailed reaction kinetic parameters [30]

The Geometrical Model and Boundary Conditions
The numerical analyses included 2 geometric configurations: A single spherical OC particle within a uniform flow field, and a random packed bed with 597 OC spheres in a cylindrical domain. For the former configuration, meshes were carefully generated in the fluid domain and the porous domain inside the sphere, which has a diameter of 5 mm in most situations and can be 3-7 mm when studying the influence of particle sizes. The default operation parameter settings are listed in Table 2, and the specified variable parameters are listed in Table 3.  For the packed bed configuration, 597 spheres with a diameter of 5 mm were randomly packed in a cylindrical column with a diameter of 30 mm (aspect ratio D/d: 6), as shown in Figure 1, and the packing model was generated with an algorithm proposed by Partopour and Dixon [34]. The height of the packing zone was approximately x = 100 mm, with an average bed voidage of 0.46. Moreover, both upstream and downstream ends of the packing zone were provided with sufficient length in order to eliminate the boundary effects. Both the inlet and outlet of the cylindrical domain were imposed Dirichlet boundary conditions, i.e., velocity-inlet and pressure-outlet in terms of Ansys Fluent terminology. On the cylinder wall, the nonslip boundary condition was set for the momentum equation, and constant temperature was set for the energy equation. Each sphere was set as a porous region, and all the conjugate fluxes were automatically ensured by the flow solver. The boundary condition settings for the first geometrical configuration of one single sphere were exactly same.
particle within a uniform flow field, and a random packed bed with 597 OC spheres cylindrical domain. For the former configuration, meshes were carefully generated in fluid domain and the porous domain inside the sphere, which has a diameter of 5 mm most situations and can be 3-7 mm when studying the influence of particle sizes. default operation parameter settings are listed in Table 2, and the specified vari parameters are listed in Table 3. For the packed bed configuration, 597 spheres with a diameter of 5 mm w randomly packed in a cylindrical column with a diameter of 30 mm (aspect ratio D/d as shown in Figure 1, and the packing model was generated with an algorithm propo by Partopour and Dixon [34]. The height of the packing zone was approximately x = mm, with an average bed voidage of 0.46. Moreover, both upstream and downstr ends of the packing zone were provided with sufficient length in order to eliminate boundary effects. Both the inlet and outlet of the cylindrical domain were impo Dirichlet boundary conditions, i.e., velocity-inlet and pressure-outlet in terms of An Fluent terminology. On the cylinder wall, the nonslip boundary condition was set for momentum equation, and constant temperature was set for the energy equation. E sphere was set as a porous region, and all the conjugate fluxes were automatically ensu by the flow solver. The boundary condition settings for the first geometrical configura of one single sphere were exactly same. The polyhedral grid was chosen for meshing with the extremely complex struc of the packed bed reactor, which has a good orthogonality and can efficiently comp the simulation works at a reasonable computational cost. In addition, an artificial ga 1.5% of the particle diameter was treated in contact point for preventing cell skewness accelerating convergence. As set forth by Eppinger et al. [16], several base sizes w The polyhedral grid was chosen for meshing with the extremely complex structure of the packed bed reactor, which has a good orthogonality and can efficiently complete the simulation works at a reasonable computational cost. In addition, an artificial gap of 1.5% of the particle diameter was treated in contact point for preventing cell skewness and accelerating convergence. As set forth by Eppinger et al. [16], several base sizes were suggested, and all mesh parameters were set to a percentage of the base sizes, thus simplifying the steps of grid independence verification. The mesh generation parameters in the current works are shown in Table 4, and the grid verification results are discussed in Section 3.2.

Results and Discussion
In this work, the heterogeneous reaction of methane CLC and its relevant conjugate transport processes inside a relatively large spherical porous particle were numerically solved based on the simplified noncatalytic reaction model [33] in order to investigate how convection in the surrounding flow field and the particle properties influence the species transport and heterogeneous reaction inside the OC particle. The effects of the Reynolds number, OC content load, internal pore size and particle size on the internal diffusion and heterogeneous reaction were investigated thoroughly. Based on the comprehensive analyses on single OC particles, the investigations were extended to the reaction and its relevant conjugate transports in a lab-scale random packed bed reactor to understand, at a typical external convection condition, how internal diffusion limitations interact with the complicated interstitial flow in the packed bed and finally determine the reactor-scale CLC processes.

Validation of Numerical Methods
The CLC that occurs in a large porous OC particle includes processes like species diffusion in the pores, gas-solid non-catalytic reactions and gas fluxes across the particle surface interacting with external convection. In the present work, those processes were numerically solved with a flow solver combined with the simplified heterogeneous reaction model which determines the species and energy sources in the flow solver.
At first, the flow solver was validated with the well-recognized correlations of convective heat and mass transfer around a sphere. The average Nusselt number around the spherical particle was compared in a range of Reynolds numbers of 1-1000 with the correlation by Whitaker [35], as in Equation (15). In the same way, the conjugate mass transfer around the particle was validated by comparing the average Sherwood number in the same range with the correlation by Bird et al. [36], as in Equation (16). Figure 2a,b clearly shows that the numerical solution of the present work is in reasonable agreement with the corrections. In addition, the numerical results of the Nusselt and Sherwood numbers were within ±4.5% and ±10% when compared to the well-known correlations, respectively.

Effects of the External Convection
It can be expected that the gas diffusion resistance must be larger than usual in a relatively large porous particle and that the rate of heterogeneous reaction inside can be limited by the relatively slow internal diffusion. On the other hand, the internal diffusion certainly interacts with the external convection in the surrounding flow field. With the

Effects of the External Convection
It can be expected that the gas diffusion resistance must be larger than usual in a relatively large porous particle and that the rate of heterogeneous reaction inside can be limited by the relatively slow internal diffusion. On the other hand, the internal diffusion certainly interacts with the external convection in the surrounding flow field. With the help of the CFD-based numerical approaches mentioned in Section 2, the influences of the external convection on the CLC gas-solid reactions and the relevant conjugate transports were quantitatively investigated.
In nature, all processes involved in CLC, including reactions and transports, are transient. Figure 3a illustrates the conversion of Cu developing during the oxidation reaction for various Reynolds numbers over time. The OC conversion evolution is obviously determined by the rate of heterogeneous reaction inside the porous particle, which, in turn, is influenced by the local concentration of gases diffused inward and outward in the porous particle. The oxygen that reacted with the Cu distributed in the particle was supplied by convection in the flowing air around the particle. According to the curves in Figure 3a, the external convection generally enhanced the internal diffusion and consequently enhanced the gas-solid CLC reactions. It can be found that the reaction was obviously slower when Reynolds number was small (Re p = 1 and 10), which implies that the external convection does not supply enough reacting species into the OC particle, that the reaction process is limited by the convection and that an increasing Reynolds number accelerates the reaction. However, in a range of higher Reynolds numbers, when the Reynolds number was higher than 50, the conversion-time curves in Figure 3a were extremely close to each other. This indicates that the reaction rate is almost same at Re p = 50, 100 and 1000, and in terms of the current particle properties, there is a threshold for influence of the external convection and the limiting factor for the gas-solid reactions becomes the internal diffusion. Figure 3b shows development of the average ∆T in the OC particle with the reaction time at various Reynolds numbers. In exothermic CLC oxidation reactions, the ∆T is determined by the heat release of reaction and the conjugate heat transfer in and around the particle. The ∆T maxima on the curves imply a balance of heat release of reaction and heat removal of external convection. These time moments also indicate how fast the reactions occur under various Reynolds numbers. Obviously, higher Reynolds numbers (Re p = 50, 100, 1000) lead to much faster reactions. On the other hand, stronger external convections also ensure much better cooling that results in obviously lower ∆T.

Effects of OC Particle Properties
Besides the external convection around an OC particle, particle properties also play a critical role in determining the heterogeneous reaction rate and species transport in a large OC particle. Therefore, it is particularly necessary to study the effects of various properties of oxygen carrier particles, such as internal pore size, oxygen carrier load and particle size, on the CLC reaction and internal diffusion.

Effects of OC Particle Properties
Besides the external convection around an OC particle, particle properties also play a critical role in determining the heterogeneous reaction rate and species transport in a large Processes 2021, 9, 125 8 of 21 OC particle. Therefore, it is particularly necessary to study the effects of various properties of oxygen carrier particles, such as internal pore size, oxygen carrier load and particle size, on the CLC reaction and internal diffusion.
Internal Pore Size of OC Particles In general, both molecular diffusion and Knudsen diffusion contribute to the effective diffusion in porous media. For porous catalytic particles with small pore sizes, Knudsen diffusion may be dominant, since the mean free path of gas molecules may be large with respect to the pore size or the microchannel radius inside the particle. For porous OC particles, it is reasonable to assume that the pore sizes are much smaller than the mean free path of the gas molecules under the conditions of CLC processes [31] and that the diffusion resistance is mainly dominated by Knudsen diffusion [30]. According to Equation (9), Knudsen diffusivity is proportional to the internal pore size and it can be expected that the internal pore size of the porous OC particles will significantly affect Knudsen diffusion and, consequently, the effective internal diffusion.
As expected, the conversion curves in Figure 4a illustrate the accelerated reaction rates by increasing the internal pore size due to the increase of the Knudsen diffusivity. Obviously, it can be found that the influence of pore size on reaction rate was not linear in the current range. For relatively large pore sizes, such as 15 mm and 20 nm, the oxidation reaction performed almost 100% conversion in 200 s, although the 20 nm case showed faster conversion of the OC content. Cases of 10 mm and 5 nm pore sizes presented much larger deviation and both reached lower conversions than the other two cases. Especially in the case of the 5 nm pore size, only 80% conversion was performed in 200 s due to its poor internal diffusion. Therefore, it can be concluded that pore size is an essential factor when choosing or designing OC supporting materials which clearly determines the internal diffusion limitation. Figure 4b illustrates the profiles of ∆T at various pore sizes over time, in which ∆T was found to increase with the pore size due to heat release of the gas-solid reaction. Since the Reynolds number of the external convection was set to be the same (meaning that the cooling flux was nearly same), the difference of ∆T for the different pore size was mainly produced by the difference of heat release. Therefore, the larger pore size of the OC particles led to more heat release, which corresponded to a higher reaction rate. It can also be seen that the curve of 5 nm pore descended obviously slower than the other three cases because of its much slower conversion.

Oxygen Carrier Load
In CLC processes, the oxygen carrier (Cu/CuO in this work) loaded in the porous inert supporter (Al2O3 in this work) directly reacts with fuel gas or oxygen in air, releasing or regenerating lattice oxygen to complete the chemical cycle. Therefore, the OC content load or concentration directly influences the CLC reactions. Three typical active OC content loads were selected here to study the influence of oxygen carrier loads on heterogeneous reaction of CLC.
In Figure 5a, the oxidation reaction conversion of a single particle at different active

Oxygen Carrier Load
In CLC processes, the oxygen carrier (Cu/CuO in this work) loaded in the porous inert supporter (Al 2 O 3 in this work) directly reacts with fuel gas or oxygen in air, releasing Processes 2021, 9, 125 9 of 21 or regenerating lattice oxygen to complete the chemical cycle. Therefore, the OC content load or concentration directly influences the CLC reactions. Three typical active OC content loads were selected here to study the influence of oxygen carrier loads on heterogeneous reaction of CLC.
In Figure 5a, the oxidation reaction conversion of a single particle at different active content loads is shown over time. It can be seen that it took longer time for cases with higher loads to perform a complete conversion, although the reaction rate was found to be quite close for all the three cases, in the early 30 s, while the overall conversion of particles with high active loads was lower at the same time. This again seems to prove that the reaction was limited by the internal diffusion, because the same reacting species flux determined by the internal diffusion consumed the OC at the same rate and therefore spent longer to finish the full conversion for the higher load cases.
(a) (b) Figure 4. The profiles of (a) overall conversion and (b) ∆T at various particle pores size over time.

Oxygen Carrier Load
In CLC processes, the oxygen carrier (Cu/CuO in this work) loaded in the porous inert supporter (Al2O3 in this work) directly reacts with fuel gas or oxygen in air, releasing or regenerating lattice oxygen to complete the chemical cycle. Therefore, the OC content load or concentration directly influences the CLC reactions. Three typical active OC content loads were selected here to study the influence of oxygen carrier loads on heterogeneous reaction of CLC.
In Figure 5a, the oxidation reaction conversion of a single particle at different active content loads is shown over time. It can be seen that it took longer time for cases with higher loads to perform a complete conversion, although the reaction rate was found to be quite close for all the three cases, in the early 30 s, while the overall conversion of particles with high active loads was lower at the same time. This again seems to prove that the reaction was limited by the internal diffusion, because the same reacting species flux determined by the internal diffusion consumed the OC at the same rate and therefore spent longer to finish the full conversion for the higher load cases. Although the oxygen carrier particles with higher OC loads were more obviously limited by the internal diffusion and the gas-solid reaction took longer to complete, an increasing OC load still means that a more intensive exothermic reaction took place inside Although the oxygen carrier particles with higher OC loads were more obviously limited by the internal diffusion and the gas-solid reaction took longer to complete, an increasing OC load still means that a more intensive exothermic reaction took place inside the particle, where more reaction heat was generated and higher ∆T was formed, as shown in Figure 5b.

OC Particle Sizes
In order to accurately understand the effect of particle diameter on the internal diffusion of large particles, several single particles with different sizes were simulated, in which the equal load density and the equal load content of active ingredient of oxygen carrier were, respectively, considered.
In the first group, the load density of single OC particles was ensured to be the same, that is, the total amount of active substances per unit volume was the same. It can be seen from Figure 6a that full conversion was quickly performed for cases of small particle sizes (d = 3 mm and 4 mm), and for the other three cases of larger particles sizes, only about 80% to 90% of conversion was achieved in the same duration. Under the same load density, it is obvious that the total amount of reaction moles needed to complete for cases of the larger the particle sizes was greater than those of the smaller ones. In addition, a reducing particle size means a decrease in the diffusion distance. Consequently, the reduced species diffusion resistance speeds up the reaction. The ∆T profile is shown in Figure 6b, although the reaction rate of the small-diameter oxygen carrier particles was accelerated due to the smaller internal diffusion resistance, which reached the maximum temperature more rapidly. The overall reaction ∆T was also lower due to the low total active content. In fact, the higher maxima of the curves for larger particle sizes also resulted from the larger resistance of heat removal due to their larger particle diameters.
80% to 90% of conversion was achieved in the same duration. Under the same load density, it is obvious that the total amount of reaction moles needed to complete for cases of the larger the particle sizes was greater than those of the smaller ones. In addition, a reducing particle size means a decrease in the diffusion distance. Consequently, the reduced species diffusion resistance speeds up the reaction. The ΔT profile is shown in Figure 6b, although the reaction rate of the small-diameter oxygen carrier particles was accelerated due to the smaller internal diffusion resistance, which reached the maximum temperature more rapidly. The overall reaction ΔT was also lower due to the low total active content. In fact, the higher maxima of the curves for larger particle sizes also resulted from the larger resistance of heat removal due to their larger particle diameters.
(a) (b) Figure 6. The profiles of (a) overall conversion and (b) ∆T with various particle size in the equal active load density over time.
In the second group, the total load content of single OC particle was ensured to be the same, which implies that the smallest particle, with a size of 3 mm in diameter, was provided with the highest load density. The temporal profiles of the conversion and ΔT for cases of different particle sizes with the same total active content are shown in Figure  7. As can be seen from the conversion profile in Figure 7a, the overall conversion of the smaller oxygen carrier particle was lower than that of the large diameter cases, and the net reaction rate was also the lower (indicated by the curve slope). It is obvious that the outer convective-diffusion area of the smaller particles was smaller, which limited the inward and outward fluxes of species diffusion and resulted in a small overall conversion and reaction rate. Similarly, due to the smaller convective area of the small particles and the larger load density, the reaction heat released in the smaller particles could not be removed in time, resulting in a sharp increase and higher peak value in the ΔT, as shown In the second group, the total load content of single OC particle was ensured to be the same, which implies that the smallest particle, with a size of 3 mm in diameter, was provided with the highest load density. The temporal profiles of the conversion and ∆T for cases of different particle sizes with the same total active content are shown in Figure 7. As can be seen from the conversion profile in Figure 7a, the overall conversion of the smaller oxygen carrier particle was lower than that of the large diameter cases, and the net reaction rate was also the lower (indicated by the curve slope). It is obvious that the outer convective-diffusion area of the smaller particles was smaller, which limited the inward and outward fluxes of species diffusion and resulted in a small overall conversion and reaction rate. Similarly, due to the smaller convective area of the small particles and the larger load density, the reaction heat released in the smaller particles could not be removed in time, resulting in a sharp increase and higher peak value in the ∆T, as shown in Figure 7b, which may have led to local hot spots and, consequently, the deactivation of the oxygen carrier in the reactor. in Figure 7b, which may have led to local hot spots and, consequently, the deactivation of the oxygen carrier in the reactor.
(a) (b) Figure 7. The profiles of (a) overall conversion and (b) ∆T with various particle sizes in the equal active load content over time.

The Reaction and Transport Characteristics in a Random Packed Bed
In Section 3.1, the reaction and transport characteristics of a single large OC particle are thoroughly investigated. The investigations are extended in this section to the reaction and its relevant conjugate transports in a lab-scale random packed bed reactor to under-

The Reaction and Transport Characteristics in a Random Packed Bed
In Section 3.1, the reaction and transport characteristics of a single large OC particle are thoroughly investigated. The investigations are extended in this section to the reaction and its relevant conjugate transports in a lab-scale random packed bed reactor to understand, at a typical external convection condition, how the intraparticle diffusion interacts with the complicated interstitial flow in the packed bed and finally affects the reactor-scale CLC processes.
The heterogeneous reactions, along with the intraparticle diffusion and interstitial transports involved in the CLC process in a fixed bed randomly packed with 597 porous spheres, were numerically simulated at Re p = 66, which is considered at or above the threshold of external convection enhancement (refer to Figure 3a), with the CFD method described in Section 2. The geometric dimensions, physical properties and operating conditions with modelling are shown in Table 5. The radial porosity distribution profile of the current packed bed is shown in Figure 8 and compared with the available references. The profile clearly shows that the present packing model is in good agreement with the correlation proposed by de Klerk [37] and the algorithm by Muller [38]. In addition, the deviations of the radial porosity distribution in the current random packing model were within ±8% (except near the axis) and ±4% for the correction by de Klerk [37] and algorithm by Mueller [38], respectively. The radial porosity distribution profile of the current packed bed is shown in Figure 8 and compared with the available references. The profile clearly shows that the present packing model is in good agreement with the correlation proposed by de Klerk [37] and the algorithm by Muller [38]. In addition, the deviations of the radial porosity distribution in the current random packing model were within ±8% (except near the axis) and ±4% for the correction by de Klerk [37] and algorithm by Mueller [38], respectively.  [37] and an algorithm from Mueller [38].
As seen in Figure 9, the profile of axial pressure drop distribution indicates that the current packing model is in good agreement with the well-known correlation proposed by Ergun [39] in Equation (17). Figure 8. Comparison of radial porosity distribution profile for the current packing model with a correction from de Klerk [37] and an algorithm from Mueller [38].
As seen in Figure 9, the profile of axial pressure drop distribution indicates that the current packing model is in good agreement with the well-known correlation proposed by Ergun [39] in Equation (17). Figure 8. Comparison of radial porosity distribution profile for the current packing model with a correction from de Klerk [37] and an algorithm from Mueller [38].
As seen in Figure 9, the profile of axial pressure drop distribution indicates that the current packing model is in good agreement with the well-known correlation proposed by Ergun [39] in Equation (17). Figure 9. Comparison of the axial pressure drop for simulation with the Ergun [39] correction at various Reynolds numbers. Figure 9. Comparison of the axial pressure drop for simulation with the Ergun [39] correction at various Reynolds numbers.
The grid independence test was performed with the single OC particle configuration. The base size of the grid was set as 3 mm, 4 mm and 5 mm, respectively, and all the validation simulations were carried out at Re p = 100. The profiles of ∆T and conversion distribution were set on the reference line at t = 50 s, as shown in the Figure 10a The grid independence test was performed with the single OC particle configuration. The base size of the grid was set as 3 mm, 4 mm and 5 mm, respectively, and all the validation simulations were carried out at Rep = 100. The profiles of ΔT and conversion distribution were set on the reference line at t = 50 s, as shown in the Figure 10a,b. Obviously, the grid independence tests worked well for all three basic sizes. Thus, all the simulations were done with a base size of 5 mm.
(a) (b) Figure 10. Grid independence test at various base sizes with a single OC particle: The profile of (a) conversion distribution and (b) reaction ∆T distribution on the reference line at t = 50 s. Figure 11 shows the distribution of velocity (left) and ∆T (right) in an axial plane at t = 150 s. From the velocity distribution, it can be seen that there was obvious channelling flow near the wall and core region, which is typical in packed beds with a small diameter ratio (D/d = 6), corresponding to the radial porosity distribution in Figure 8. In addition, there was a great temperature gradient in radial direction, which was caused by the constant cool wall boundary (T = 873 K) and the enhanced convective heat transfer by the channelling flow near the wall or the so-called wall effects.  flow near the wall and core region, which is typical in packed beds with a small diameter ratio (D/d = 6), corresponding to the radial porosity distribution in Figure 8. In addition, there was a great temperature gradient in radial direction, which was caused by the constant cool wall boundary (T = 873 K) and the enhanced convective heat transfer by the channelling flow near the wall or the so-called wall effects.
(a) (b) Figure 10. Grid independence test at various base sizes with a single OC particle: The profile of (a) conversion distribution and (b) reaction ∆T distribution on the reference line at t = 50 s. Figure 11 shows the distribution of velocity (left) and ∆T (right) in an axial plane at t = 150 s. From the velocity distribution, it can be seen that there was obvious channelling flow near the wall and core region, which is typical in packed beds with a small diameter ratio (D/d = 6), corresponding to the radial porosity distribution in Figure 8. In addition, there was a great temperature gradient in radial direction, which was caused by the constant cool wall boundary (T = 873 K) and the enhanced convective heat transfer by the channelling flow near the wall or the so-called wall effects. Figure 11. The distribution of velocity (left) and temperature (right) in an axial plane at t = 150 s. Figure 12 illustrates the distribution of O2 volume fraction, CuO mole fraction and temperature inside particle in an axial plane at various time moments. In Figure 12a, the Figure 11. The distribution of velocity (left) and temperature (right) in an axial plane at t = 150 s. Figure 12 illustrates the distribution of O 2 volume fraction, CuO mole fraction and temperature inside particle in an axial plane at various time moments. In Figure 12a, the O 2 volume fraction distribution is shown at several moments, which presents the evolution of the O 2 concentration field in the bed. It can be found that the O 2 concentration gradient, larger or smaller, existed inside the particles for all the time moments except t = 300 s. The local species concentration in particles depends on both the inward diffusion flux and the consumption rate due to the local reaction. The latter is determined by both the local species concentration and the local solid OC content, according to the Arrhenius equation. Therefore, although a complicated mechanism determines the local O 2 concentration inside a particle, the influence of internal diffusion resistance is essential or even dominant to the local species concentration, which restricts oxygen access to the particle core and thereby limits the local reaction rate in the particle centre. The most obvious presentation of the internal diffusion limitation was the concentration gradient in regions close to the bed inlet and outlet in the earlier times (20 s, 50 s and 100 s), which was characterized by the "orange rings" in the upstream (lower part of the column) and the "blue spots" in the downstream (upper part). The higher concentration layers near the particle surfaces (the "orange rings") in the inlet region were built up due to the fast depletion of Cu and easy O 2 diffusion near the particle surface, which could not easily occur further inside the particles during the earlier period. In the downstream region, O 2 was not sufficiently provided due to the consumption in the upstream. O 2 slowly diffused and reacted with Cu in the particles and was finally used up before reaching the particle cores, leaving the cores "blue". The "orange ring" type of high gradient distribution evolved along the axial direction with time and existed even at t = 200 s. To fully demonstrate the effects of the bed wall and the radial porosity distri the scalar distributions of O2 volume fraction, CuO mole fraction and ΔT in a radiall plane (x = 50 mm in middle of packing zone) are shown at various times in Figure 1 small aspect ratio (D/d) packed bed, the current reactor showed obvious wall eff illustrated in Figure 11. That means there was stronger convection near the wall the core region, which is proved by the higher near-surface concentration of O2 an of the outer particles radially in Figure 13a,b. Although the condition of Rep = 66 The distribution of the mole fraction of CuO can be regarded as the local conversion distribution, which is closely related to the oxygen concentration distribution, as shown in Figure 12b. When compared with the O 2 distribution in Figure 12a, the CuO distribution is quite similar but does not coincide completely with the former. In general, CuO formed gradually from upstream to downstream and from the particle surface to the interior. The low O 2 concentration in the particle cores caused by the internal diffusion limitation resulted in a small reaction rate, and thereby reduced the local conversion. As time went on, the Cu near the surface of particle was gradually consumed, and O 2 reached the particle core, where it reacted with Cu. It can be found that the complete conversion of the oxygen carrier Cu in packed bed took a long time, and even at t = 200 s, unreacted oxygen carrier remained in the cores of the particles downstream.
As heat generation is an essential goal of CLC, the investigation of ∆T resulting from the reaction is of great significance to the design and development of packed bed reactors. In Figure 12c, the ∆T distribution is shown at several moments. Basically, the ∆T in th packed bed reactor was affected by the rate of exothermic reaction and convective heat removal. Furthermore, the ∆T in the packed bed was also determined by the effect of temperature superposition caused by heat flow from upstream. At t = 20 s and 50 s, the temperature increased with time, and at t = 100 s, a hot zone began to form due to the effect of the heat flow from upstream. Over time, the hot zone moved along the flow direction, and at t = 300 s, most parts of the hot zone moved out of the bed. After this moment, the packed bed started to substantially cool from the upstream.
To fully demonstrate the effects of the bed wall and the radial porosity distribution, the scalar distributions of O 2 volume fraction, CuO mole fraction and ∆T in a radially sliced plane (x = 50 mm in middle of packing zone) are shown at various times in Figure 13. As a small aspect ratio (D/d) packed bed, the current reactor showed obvious wall effects, as illustrated in Figure 11. That means there was stronger convection near the wall than in the core region, which is proved by the higher near-surface concentration of O 2 and CuO of the outer particles radially in Figure 13a,b. Although the condition of Re p = 66 would imply a diffusion-controlled process if a single particle was considered (as discussed in Section 2), it seems that the stronger convection near the wall still enhanced the across-surface species transport inward and the reaction in the layer near the particle surface. On the other hand, the local conversion in the core region of the particles was very low for both inner and outer particles due to the internal diffusion limitation. In Figure 13c, it is clearly observed that the radial temperature gradient developed from the beginning to the end of the process, which was caused by both the cool wall and wall effects. Figure 14 shows the overall conversion profile in the packed bed over time. Although there was an axial delay of the reaction and species transport in the packed bed, the trend of the conversion-time curve in Figure 14 is quite similar to those of the higher Reynolds numbers in Figure 3a. The curve was somehow linear in the first half of the whole process, but the conversion became obviously slower after 150 s. Obviously, due to the large diffusion area and short diffusion distance in the layer near the particle surface, the internal diffusion capacity was much higher than that in the particle core. Therefore, the reaction rate was obviously higher during 0-150 s than later, which explains why the conversion during the second half of the process came slower. By more careful comparison, it can be found that the full conversion for the packed bed was almost 100 s longer than that in the case of the single particle. The two cases are certainly different considering the extremely complex interstitial channels in the packed bed. Except the difference of the local flow fields between the two cases, the shortage of O 2 downstream in the bed due to the consumption of O 2 can result in a smaller diffusive and convective flux into the particle cores downstream, which can determine a slower conversion.
imply a diffusion-controlled process if a single particle was considered (as discussed in Section 2), it seems that the stronger convection near the wall still enhanced the across-surface species transport inward and the reaction in the layer near the particle surface. On the other hand, the local conversion in the core region of the particles was very low for both inner and outer particles due to the internal diffusion limitation. In Figure 13c, it is clearly observed that the radial temperature gradient developed from the beginning to the end of the process, which was caused by both the cool wall and wall effects.  In order to quantitatively illustrate how the temperature field evolves along the axial direction, the curves of ∆T over time at various axial positions are shown in Figure 15. The oxidation process in the packed bed that we simulated was an exothermic reaction, in which each axial position underwent a heating from the reaction heat and convective heat flux from its upstream, and a cooling flux contributed by the cool bed wall. The complex superposition of the heat fluxes determines different history curve of ∆T at various axial positions. Among the curves, the one at x = 45 mm showed up the maximum peak value, which was located approximately at the middle of the bed. Although the maxima of x = 75 and 95 mm were somewhat lower, their relevant curves also maintained higher values during the whole process. Obviously, the downstream length starting from the middle of the bed received sufficient heat flux from its upstream, so that the ∆T quickly reached and maintained a higher value. However, the reaction that occurred in that part could have been delayed and slow. By contrast, the maximum temperature close to the reactor inlet (0-25 mm) reached the peak value much faster but the maxima were substantially smaller than those at positions behind them. Obviously, at the former positions, the heat flux from their upstream was not yet sufficiently large. In fact, the value of maximum ∆T at the axial positions mainly reflects how much the convective heat flux from upstream was superposed on the local heat generation. The latter is almost the same, although there was a delay in the downstream direction. Figure 16 illustrates the local maximum ∆T (left) and the time moment to reach the maximum (right) in the axial positions, which quantitatively help us understand the distribution of hot spots in axial positions over time. Obviously, the local maximum ∆T for almost 80% length of the bed (x = 20-100 mm) was beyond 100 K, and this maximum was reached in around 200 s.
Processes 2021, 9, x FOR PEER REVIEW 17 of 22 Figure 14 shows the overall conversion profile in the packed bed over time. Although there was an axial delay of the reaction and species transport in the packed bed, the trend of the conversion-time curve in Figure 14 is quite similar to those of the higher Reynolds numbers in Figure 3a. The curve was somehow linear in the first half of the whole process, but the conversion became obviously slower after 150 s. Obviously, due to the large diffusion area and short diffusion distance in the layer near the particle surface, the internal diffusion capacity was much higher than that in the particle core. Therefore, the reaction rate was obviously higher during 0-150 s than later, which explains why the conversion during the second half of the process came slower. By more careful comparison, it can be found that the full conversion for the packed bed was almost 100 s longer than that in the case of the single particle. The two cases are certainly different considering the extremely complex interstitial channels in the packed bed. Except the difference of the local flow fields between the two cases, the shortage of O2 downstream in the bed due to the consumption of O2 can result in a smaller diffusive and convective flux into the particle cores downstream, which can determine a slower conversion. In order to quantitatively illustrate how the temperature field evolves along the axial direction, the curves of ∆T over time at various axial positions are shown in Figure 15. The oxidation process in the packed bed that we simulated was an exothermic reaction, in which each axial position underwent a heating from the reaction heat and convective heat flux from its upstream, and a cooling flux contributed by the cool bed wall. The complex superposition of the heat fluxes determines different history curve of ∆T at various axial positions. Among the curves, the one at x = 45 mm showed up the maximum peak value, which was located approximately at the middle of the bed. Although the maxima of x = 75 and 95 mm were somewhat lower, their relevant curves also maintained higher values during the whole process. Obviously, the downstream length starting from the middle of the bed received sufficient heat flux from its upstream, so that the ∆T quickly reached and maintained a higher value. However, the reaction that occurred in that part could have been delayed and slow. By contrast, the maximum temperature close to the reactor inlet (0-25 mm) reached the peak value much faster but the maxima were substantially smaller than those at positions behind them. Obviously, at the former positions, the heat flux from their upstream was not yet sufficiently large. In fact, the value of maximum ∆T at the axial positions mainly reflects how much the convective heat flux from upstream was superposed on the local heat generation. The latter is almost the same, although there was a delay in the downstream direction. Figure 16 illustrates the local maximum ∆T (left) and     Similarly, a temperature gradient formed and evolved in the radial direction during the oxidation process. Figure 17 shows that the temperature distribution at different radial coordinate points on the radial plane (x = 45 mm) developed with time. In fact, the radial temperature gradient was built up, mainly due to the cooling flux toward the cool wall and the axial convective flux, because of the channelling flow in the high-voidage region near the bed wall. However, the local heat generation may have been higher to some extent because of the enhanced species transport due to channelling flow.
In this paper, a typical small diameter ratio (D/d = 6) was selected to randomly generate a packed bed reactor with a diameter of 30 mm and a height of about 100 mm. The Similarly, a temperature gradient formed and evolved in the radial direction during the oxidation process. Figure 17 shows that the temperature distribution at different radial coordinate points on the radial plane (x = 45 mm) developed with time. In fact, the radial temperature gradient was built up, mainly due to the cooling flux toward the cool wall and the axial convective flux, because of the channelling flow in the high-voidage region near the bed wall. However, the local heat generation may have been higher to some extent because of the enhanced species transport due to channelling flow. numerical simulation results give us a comprehensive understanding of the temporal evolution and spatial distribution of internal scalars. In addition, changing the diameter of the reaction tube and the axial packing height in a reasonable range still conform to the axial and radial distribution of our conclusion.

Conclusions
In this work, a transient heterogeneous noncatalytic gas-solid reaction and its relevant conjugate transport processes in and around large oxygen carrier particles were numerically investigated with a discrete particle CFD method. A simplified heterogeneous non-catalytic gas-solid reaction model was integrated into the three-dimensional flow solver to reduce the computational costs of intraparticle modelling.
Numerical analyses of a large single OC particle were first carried out to understand the effects of external convection and particle properties on the intraparticle diffusion characteristics and conjugated transports in the process of CLC based CuO. The analyses indicate that there is a threshold for influence of the external convection roughly at or In this paper, a typical small diameter ratio (D/d = 6) was selected to randomly generate a packed bed reactor with a diameter of 30 mm and a height of about 100 mm. The numerical simulation results give us a comprehensive understanding of the temporal evolution and spatial distribution of internal scalars. In addition, changing the diameter of the reaction tube and the axial packing height in a reasonable range still conform to the axial and radial distribution of our conclusion.

Conclusions
In this work, a transient heterogeneous noncatalytic gas-solid reaction and its relevant conjugate transport processes in and around large oxygen carrier particles were numerically investigated with a discrete particle CFD method. A simplified heterogeneous non-catalytic gas-solid reaction model was integrated into the three-dimensional flow solver to reduce the computational costs of intraparticle modelling.
Numerical analyses of a large single OC particle were first carried out to understand the effects of external convection and particle properties on the intraparticle diffusion characteristics and conjugated transports in the process of CLC based CuO. The analyses indicate that there is a threshold for influence of the external convection roughly at or above Re p = 50, beyond which the limiting factor for the gas-solid reactions becomes the internal diffusion in terms of the particle properties that have been used. On the other hand, particle properties are obviously essential factors when designing OC particles, because they can directly influence intraparticle diffusion and further determine the OC reaction performance. It was found that relatively larger internal pore sizes help to enhance the reaction rate due to their lower Knudsen diffusion resistance. The influence of active content load proves that the reaction is mainly limited by the internal diffusion, because cases of higher load density take a longer time to finish the full conversion at the same reacting species flux determined by the internal diffusion. The influence of particle sizes depends on how the amount of active OC ingredient load is allotted based on the same density of load or the same total amount for all the particles with various sizes. In the former case, the full conversion of OC for smaller particles obviously requires shorter time, because their total amount of reaction moles needed to complete is smaller under the same load density. In the latter case, the reaction for smaller particles is obviously slower, since their smaller surface convective-diffusion area limits the inward and outward fluxes of species diffusion.
To understand how the intraparticle diffusion interacts with the complicated interstitial flow in a packed bed at a typical external convection condition, the numerical investigation was extended to the reaction and its relevant conjugate transports in a labscale random packed bed CLC reactor of 597 spheres. The transient evolution of gas species concentration, OC content and temperature rise were thoroughly analysed. During the oxidation process, CuO formed gradually from upstream to downstream and from the particle surface to the interior. It was found that the full conversion for the packed bed was almost 100 s longer than that in the case of the single particle because the local transport was obviously different due to the extremely complex interstitial channels in the packed bed. In addition, the shortage of O 2 downstream in the bed due to the upstream consumption of O 2 can result in a smaller diffusive and convective flux into the particle cores downstream, which can determine a slower conversion.  Data Availability Statement: Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no competing financial interest.