Heat Transfer Numerical Analysis of a Thermochemical Solar Reactor

Corresponding Author: E.G. Espinosa-Martínez, Área de Ingeniería y Recursos Energéticos, Universidad Autónoma MetropolitanaIztapalapa, D.F. México Email: yurihillel@gmail.com Abstract: In this study a 2D transient temperature distribution in a onedimensional thermochemical solar reactor of gas-particles is presented. The gas-particles medium is contained in a reactor where an endothermic chemical reaction of some oxide-metal at high temperature is expected to be occurred. The reactor receives the concentrated solar power on the reactor walls, or directly to the gas-particles medium. As expected the temperature distribution in the reactor is quite homogeneous when it operates to close steady state conditions, but also due to the fact that the inlet cold gas is reheated with the exit warm gas. However, when the reactor starts to operate, the temperature distribution is non-homogeneous depending on how the reactor receives the concentrated solar power.


Introduction
The urgent need of taking advantage of the renewable energy sources in order to replace fossil fuels, impels to review and optimize the systems design that are used in the concentration and transformation of the thermal solar energy. The technical and economic feasibility of a Solar System for industrial use depends basically on the thermal efficiency. The solar processes at high temperatures are one of the alternatives to carry out materials transformation for the chemical industry, like calcium carbonate calcination of lime or cement, production of activated charcoal, among others. The production of energy carriers, as the hydrogen, is a main application which is focused for hydrocarbons substitution, but with very low or null gas greenhouse emissions.
The thermochemical systems additionally include the production of a fuel like storage media of the solar energy. This system has the advantage of producing long term storable energy carriers from the solar energy. The thermochemical conversion also enables solar energy transportation from Sunbelt to remote population centers. Hydrogen is one of the suitable fuels that has null impact on the environment. Therefore, the adaptation to the particularities of solar energy to the solar thermochemical processes has important "R and D" tasks in this field; thermochemical hydrogen production through water splitting, steam or CO 2 gasification of coal, solar reforming of methane, among others, are important industrial processes for producing useful chemical fuels such as hydrogen, synthesis gas and C 2 hydrocarbons. In the last years, efforts have been focused in this direction (Abanades and Flamant, 2005;Müller and Flamant, 1998;Z'Graggen et al., 2008;Steinfeld, 2002;Kodama, 2003).
The efficient production of solar hydrogen is a challenge. The system is based on the use of concentrated solar radiation as the energy source for performing high temperature reactions that produce hydrogen. Several paths are possible: From transformation of various fossil and non-fossil fuels via different routes such as natural gas steam reforming, that produce syngas, water splitting thermochemical cycles, gasification and pyrolysis of biomass that produce syngas and then, hydrogen with separation of carbon dioxide (Funk, 2001;Steinfeld, 2005;. The concentrated solar radiation can also be stored by means of chemicals bonds. Some salts as sulfates, carbonates, hydroxides and ammonia salts are reversible reactions were thermal energy can be stored (Prengle and Sun, 1976;Wentworth and Chen, 1976;Ervin, 1977;Romero-Paredes et al., 2006).
Naturally the ideal row material to produce hydrogen is water, due to its abundance, low value and the absence of CO 2 emissions during its dissociation. Thermolysis is not favored thermodynamically. It needs high temperature to have a good rate of reaction, less than 4000 K takes it out of the range for solar applications. Nevertheless the current research in the world for hydrogen production is focused on the so called redox pair cycles for water splitting (Agrafiotis et al., 2007;Steinfeld, 2005;Steinfeld et al., 1999;Ambriz et al., 1982;Sibieude et al., 1982;Valdés-Parada et al., 2011). Those cycles bypass the hydrogen and oxygen separation problem and allows to operate at lower temperatures.
Solar thermochemical engineering has emerged like an activity from multidisciplinary investigation in the laboratories or national research centers, in the industry and universities anywhere in the world (Fletcher, 2001;Palumbo et al., 2004). The challenges of this activity in the field of solar chemistry include the development of advanced solar concentrators, synthesis of active materials for the solar chemical reactions (catalytic, material redox), design and operation of new receivers and reactors able to catch and transform the concentrated solar radiation with a high thermal and exergetic performance (Agrafiotis et al., 2007). The medullar center of a chemical process is the reactor which has to be reliable and must support the mechanisms of heat transfer and chemical processes with very high efficiency. Additionally the materials, the design and the method of construction of the reactor that has to work at high temperatures, above of 1200 to 2500 K, the reactor must support the thermal shocks when using the concentrated solar energy.
The solar reactors were used from the beginnings of the research in materials with solar concentrators, by the year 1964 and since then it has been proposed an ample range of designs and materials (Funk and Reinstrom, 1964;1966;Funk, 2001;Bernard et al., 1978;Ambriz and Romero-Paredes et al., 2006;Beghi, 1986). The solar reactors generally can be classified in to two categories based on the mechanism of solar radiation catching: Those of direct coupled receiving and those of indirectly coupled receiving. A type of reactor of direct coupled receiving is the designed by Haueter et al., (1999), which consists of a rotatory cavity with a transparent window of the receiving center which is aligned with the particle set that are held by centrifugal force. With this arrangement, the solid particles are directly exposed to high-flux solar irradiation. In addition it serves simultaneously like a radiant absorber, thermal insulation and chemical reactant (Steinfeld, 2005;Abanades et al., 2007). The "confined tornado" solar reactor was designed to prevent window contamination with a minimum screen flow under the window at the aperture of the receiver (Kogan et al., 2007). The vortex solar reactor is another example of a novel reactor based on the direct irradiation on the chemical reactants (Steinfeld et al., 1998a;. Fluidized bed reactor directly exposed to high-flux solar irradiation has also been studied previously (Müller and Flamant, 1998a;Müller et al., 2003;Villafán-Vidales et al., 2012). The direct irradiation of particles suspension provide efficient means of heat and mass transfer to held the reaction. Those kind of reactors are in general more efficient that those based on indirect Irradiation mode, due that they bypass the limitations imposed by indirect high temperature heat transfer (VonZedtwitz et al., 2006;Steinfeld et al., 1997;Kraupl and Steinfeld, 2003;Hirsch et al., 2004a). The "two-cavity" solar reactor based on indirect irradiation of reactants has been used for the carbothermal reduction of Z n O (Steinfeld et al., 1998a;. This type of reactor has the advantage of eliminating the need for a transparent window, but it has the disadvantage of the limitations imposed by the materials of construction of the reactor walls, that limits the maximum operating temperature, thermal conductivity, irradiance-absorptance, resistance to thermal shocks, inertness and thermal inertia. A challenge of the solar thermochemical engineering is the optimized design of a solar reactor. Then the design and operation of novel receive/reactor systems capable of collecting efficiently the concentrated solar radiation and to perform efficiently the high temperature reaction is the main objective of this field. In order to scale up solar reactor parameters as the reactor volume, maximum working temperature, physicochemical properties of the construction materials and the reactants and temperature distribution have to be optimized, taking into account the heat transfer characteristics, the reaction rates and the transient phenomena due to the random conditions of the solar flux. Relatively few studies have deal with solar reactor modeling and simulation (Abanades et al., 2007;Agrafiotis et al., 2007;Petrasch and Steinfeld, 2007;Kogan et al., 2007;Palumbo et al., 2004;Meier et al., 1996). Some of those studies have been focused on radiative heat transfer within particle suspension exposed to concentrated solar radiation, included, steady state models based on the discrete ordinates methods, Monte Carlo (MC) methods and transient models based on the MC and Rosseland approximations and diffusion approach (Evans et al., 1987;Miller and Koenigsdorff, 1991;Mischler and Steinfeld, 1995;Hirsch and Steinfeld, 2004b;Osinaga et al., 2004;Lipinski and Steinfeld, 2004;Dombrovsky, 2002;Dombrovsky et al., 2007;Klein et al., 2007;Müller and Steinfeld, 2007). Other models and simulations have been developed in order to describe the temperature profile in the reactor and to have a multiphase flow method in the case of solid particle dissociation reaction (Abanades et al., 2007) and for the description of the thermochemical cycles for hydrogen production inside a honeycomb monolithic reactor under the influx of solar radiation on its front face (Agrafiotis et al., 2007;2005). A cavity receiver with a tubular absorber using concentrated solar energy was studied for thermochemical processes, in particular for ZnO dissociation (Melchior et al., 2008).
In this study we presents simulations of the twodimensional temperature distribution in a vertical bed reactor fixed made of particulate of an oxide metal, which uses some inert gas as a working fluid. The system receives a concentrated solar power on the reactor walls, or directly to the gas-particles medium. The fixed bed reactor model considers only an incipient state and the net particle motion is just neglected with respect to the gas velocity. Additionally, the scheme presented in this study can be applied to similar systems.

Thermochemical Reactor
A thermochemical solar reactor is an option for producing hydrogen fuel with "clean" technology. The main idea behind this type of reactor is to enhance a two steps chemical reaction of an oxide metal. One step consists of an endothermic reaction for oxide dissociation at high temperature using concentrated solar power as the heat source; the other step consists of an exothermic reaction for water reduction producing hydrogen.
In this study, we propose a thermochemical vertical reactor, whose analysis consider a fixed bed particles, which uses some inert gas as working fluid. The system receives a concentrated solar power on the reactor walls, or directly to the gas-particles medium. We have considered one-dimensional gas velocity, with a negligible net particle motion with respect to the gas velocity. We have simulated the two-dimensional temperature distribution inside the reactor due to the conduction-convection modes of heat transfer. The radiative heat transfer is considered as the heat source term in the simulations.

Mathematical model
We have considered a vertical reactor of rectangular cross-section, Fig. 1. With this geometry is natural to use a Cartesian coordinate system (x, y, z). For this problem the gas flow passes through the particles in the axial direction, z.
The thermal energy equation to describe the temperature distribution in the reactor is based in a model of an average equation or thermodynamic equilibrium (Espinosa-Paredes et al., 2013a;2013b). In our approach the medium is considered isotropic and homogeneous. For simplicity the dispersion and compressibility effects are just neglected. The simplified mean equation for the temperature reads Equation 1: where, (T) is a representative temperature of the liquidsolid phases, (Φ) is the homogeneous thermal source per volume unit which is not considered when the systems receive the concentrated solar radiation through the walls, (ρCp)g is the product of the gas density, ρ and heat capacity, C p .k ef represents the effective thermal conductivity of the medium defined by Equation 2: where, k g (s) is the gas (g) or particle (s) conductivity.
The effective property (ρCp) is defined by Equation 3: The gas velocity in the heat transfer model is given by Equation 4: where, ρ g and ρ s are the gas and particle density respectively, ε g .

Boundary Conditions
The boundary conditions for temperature are established depending how the system receives the concentrated solar radiation: (1) Through the reactor walls, or (2) directly to the gas-particles medium. For the first case, the temperature at the reactor walls is just assumed as a constant Equation 5: The temperature T wall , may be interpreted as an average temperature of the solar temperature distribution of the concentrated solar "spot" on the wall.
For the second case, the reactor walls are considered as insulated, then the condition for temperature is: For this system the concentrated solar radiation arrives at the upper side of the reactor, that is, at z = 4R where the gas exits as shown in Fig. 1. In our model the solar radiation intensity is assumed to be a Gaussian-like distribution which decay exponentially along the reactor. The radiation is considered as a source term in Equation 1.
Additionally, in both systems it is assumed that the inlet gas temperature is heated up with the recovery heat from the exit gas. In this sense the inlet gas temperature is defined as an average of the exit gas temperature distribution. Concerning the initial condition, an uniform constant temperature T 0 of the systems is considered.

Dimensionless Model
We proceeded to adimensionalize Equation 1, boundary conditions 5 to 6 and the initial conditions as follow: The length scales are measured in terms of the reactor hydraulic radius R, the time scale in terms of R 2 (ρCp)/k ef ; and the temperature scale in terms of T wall , for constant solar heat flux at the wall, or in terms of R 2 q 0 /k ef , for direct solar heat flux into the medium. The quantities T wall and q 0 are the constant wall temperature and maximum intensity of the concentrated solar radiation, respectively. The last scaled quantities are used in Equation 1, in Equation 5 to 6 and in the initial condition. The resulting dimensionless equation for temperature is Equation 7: The new variables θ, Φ, τ, X and Z correspond to the dimensionless variables of temperature, heat source, time, x and z direction respectively. We have defined the Peclet number as Pe = u g R/α ef , in terms of µ g and an effective thermal diffusivity defined byα ef = k ef /(ρCp). The factor X is only the ratio between ε g (ρCp) g and (ρCp).
The dimensionless equations of the boundary and initial conditions, for constant temperature on the walls are Equation 8: Respectively, for the case of solar heat flux directly to the medium, the respective boundary and initial conditions are Equation 10:

Numerical Solution
In order to obtain the temperature distribution in the solar reactor as a function of time and position, a transient model based on finite differences in the implicit scheme is used. The choice of this discretized scheme is common for simple geometry configurations, making it easier and/or faster than other methods (e.g., method of characteristic). We applied a first-order downstream implicit scheme for partial derivatives and a first-order upstream for time derivatives. The resulting set of equations in discretized form can be written as Equation 12: Where Equation 13 to 17: , , , In these equations the superscripts τ and τ + ∆τ indicate that the dependent variables are calculated at a previous dimensionless time step and at the next one, respectively. The subscripts i and j, denote the cell number where the variable is calculated as illustrated in Fig. 2, for purpose discussion, we consider only a 10 by 10 grid system (Espinosa-Paredes et al., 2012). The variables with subscripts j-1, i-1 and superscript τ are known, since these are the inlet variables (i.e., boundary conditions) and the initial condition, respectively.
As a result of applying the equation recursively for i = 1, 2,…,n and for j = 1, 2,…,m (with n and m integer numbers), we get a set of algebraic equations, which can be written in matrix form as Equation 20: where, A is the coefficient matrix which is pent a diagonal, x is the unknown vector and b is the coefficient vector. The numerical solution of Equation 7 was obtained using an algorithm based on the factorization of a matrix using a version of the Gaussian elimination with partial pivoting.

Numerical Validation
In order to validate our numerical model we have compared analytical solutions, described in Appendix A, against numerical solutions computed from our numerical code. Basically, four solutions to Equation 7 with specific boundary and initial conditions were tested. A numerical analysis of mesh independence was also performed. These analyses are described in this section.

Diffusion Model
In the first validation, a transient heat diffusion model without heat source is considered. This implies that the convective term XPe∂θ/∂Z and source term Φ of Equation 7 are zero. In order to solve this problem, it is enough to established four spatial boundaries and one temporal boundary. Three of the spatial boundaries are fixed to some constant temperature and the missing one is considered as isolated. Initially the resting medium is at some constant temperature, evidently this temperature has to be different from temperature of the spatial boundaries if some enhanced diffusion is desired. Figure  3a "temporal photo" of the temperature distribution along the dimensionless width X at the middle of the reactor, that is for a horizontal cut at Z = 2, is shown for both solutions: Analytical and numerical. The same "temporal photo" is plotted in Fig. 3b, excepting that the temperature distribution along the dimensionless depth Z for a cut at X = 1 (corresponding to the symmetry line in this direction) is shown. In both pictures a very good agreement between both solutions for a 21×21 mesh is observed. In fact, the maximal numerical error is around 2%.

Diffusion-Convection Model
The second validation considers a conductiveconvective transient heat model without heat source. This implies that in Equation 7 the term given by Φ is just set to zero. We have considered the same spatial and temporal boundaries as in previous test. Figure 4, as previously, we have also plotted two "temporal photos" of the temperature distribution along the dimensionless width X for a horizontal cut at Z = 2 and for the temperature distribution along the dimensionless depth Z for a cut at X = 1. In these figures a good agreement between both solutions for a 21×21 mesh is observed, however, the maximal numerical error is around 4%, twice larger than the error in the first test.

Diffusion Model with Heat Source
We have considered in this test, the same time dependent diffusion problem as in section 3.1, except that a given heat source term Φ is now included in the model. As in previous cases, we have fixed three spatial boundaries to some constant temperature and stated one spatial boundary as isolated. Similarly, an initial temperature condition of the resting medium is considered. Figure 5, we also show two "temporal graphs" of the temperature distribution: One along the dimensionless width X for a horizontal cut at Z = 2 and the other along the dimensionless depth Z for a cut at X = 1. In these figures a very good agreement between both solutions for a 21×21 mesh is observed. In fact, this is the best fitting comparison, where the numerical error is less than 1%.

Diffusion-Convection Model with Heat Source
As a final comparison, we have considered the whole problem, that is, a conductive-convective transient heat model with heat source given in Equation 7. The boundary and initial conditions remain the same as in previous cases. Figure 6, similar "temporal photos" of the temperature distribution are shown. Similarly, a good agreement between both solutions can be observed and the maximal numerical error obtained is around 4%.

Mesh Analysis
In this part a mesh independence analysis was performed. We have compared the temperature profile obtained from three different mesh with respect to the analytical solution. The test model corresponds to the diffusion-convection problem without heat source. Figure  7, the temperature profile in X direction for a horizontal cut at Z = 2 and a given time is shown. From the figure we can state that numerical solutions have a good approach to the analytical one; however, a better accuracy is obtained when the number of nodes is increased.

Results and Discussions
We have obtained numerical results of the model given by Equation 7 to 9 for a range of values of XPe running from 0 to 5 and for a range of values τ running from 0 (earlier times) to 2 (long times). Two cases were considered: Constant temperature at the wall and concentrated solar radiation directly to the medium.  Figure 8a only heat diffusion is considered while in Fig. 8b heat convection is added to the diffusion system. In both pictures, we can observe that at earlier times the temperature distribution is not homogeneous in the X direction, but just at long times when the wall temperature is closely reached. If we compare both systems, the convection-diffusion system reaches more rapidly a homogeneous temperatures distribution for the same dimensionless time values.

Constant Temperature at the Wall
In Fig. 9, a "temporal graph" of dimensionless temperature profiles along the X direction are shown. Each curve corresponds to different gas flow rates, that is, for different values of XPe. The idea behind this picture is to show the convection effect on the temperature magnitude and distribution in the reactor. We can observe that when gas flow rate is increased, higher values of temperature are more rapidly reached. For this system a maximum difference of 5% is found between the minimum value of temperatures at XPe = 0 and XPe = 5.
In Fig. 10, contour plots of temperature distribution for a fixed dimensionless time (τ = 2 are shown. In Fig.  10a corresponds to the diffusion system with fixed "cold" temperature at the gas inlet (z = 0). This plot shows homogeneous temperature distribution (close to that at the wall) in the upper half of the reactor, while in the lower half part, a non homogeneous distribution is observed due to strong diffusion effects. For the diffusion-convection case (Fig. 10b) a global homogeneous distribution of temperature is obtained along the reactor. This due to the convection effect but also to the fact that the inlet temperature of the gas interchange heat with the exit gas, increasing the inlet temperature with time. Then for "long times", the convection diffusion system achieved an homogenized temperature value close to that at the wall. Figure 11 curves of the temperature distribution along the dimensionless width (X) for the case of concentrated solar flux directly to the gas-particles medium for different values of dimensionless time are represented. The shown picture view is from an horizontal plane at Z = 2, that is, at the middle of the reactor. Figure 11a only heat diffusion is considered while in Fig. 11b a heat convection-diffusion system is analyzed. In both pictures, we can observe that at earlier times (τ close to zero) the temperature distribution is less homogeneous than the temperature distribution at long times (τ close to two). However, if we compare both systems, the convection-diffusion system reaches more rapidly an homogeneous temperatures distribution for same dimensionless time values. In any case temperatures distribution in this system is practically more homogeneous with respect to the fixed wall temperature system. Figure 12, a "temporal graph" of dimensionless temperature profiles along the X direction are shown. Each curve corresponds to different gas flow rates. As the previous system, the convection effect on temperature intensity and distribution inside the reactor is shown. As expected temperature values increase when the gas flow rate is increased. Also for this system a maximum difference of 5% is found between the minimum value of temperatures at XPe = 0 and XPe = 5. Figure 13, contour plots of temperature distribution for a fixed dimensionless time (τ = 2) are shown. Figure  13a corresponds to the diffusion system with fixed "cold" temperature at the gas inlet (z = 0). This plot shows practically homogeneous temperatures distribution along the direction, but not along the direction where a high temperature gradient exists. For the diffusion-convection case (Fig. 13b) an homogeneous temperature distribution is obtained at the lower quarter of the reactor (between Z = 0 and Z = 0). This means that the heated inlet gas, by the gas outflow, is mostly convected in this region. On the other hand, at the upper quarter of the reactor (between Z = 3 and Z = 4) a homogeneous temperature distribution exists. This homogeneity in temperature is due to the fact that the maximum solar heat source concentration is at the center (at X = 1) and the convective effects from the inlet has not raised to that region. The visible dark blue area in Fig. 13b is a memory effect of the temperature at previous times.

Summary and Conclusions
In this study we have developed a model for describing the two-dimensional temperature distribution in a vertical bed reactor fixed made of particulate of an oxide metal, which uses some inert gas as a working fluid. The system receives a concentrated solar power on the reactor walls, or directly to the gasparticles medium. The fixed bed reactor model considers only an incipient state and the net particle motion is just neglected with respect to the gas velocity. The theoretical model for temperature is based on an average equation or thermodynamic equilibrium. The medium is contained in a reactor where an endothermic reaction of oxide dissociation at height temperature is expected to be occurred.
This model could be interpreted as gas flow through a porous media, where the porosity of the medium depends of the flow itself. Two types of reactor configuration were studied depending how the concentrated solar power heats the system: (1) Through the reactor walls and (2) directly to the gasparticles medium. In our model, radiation heat transfer is considered only as the given heat source term in the equations. Then, the radiation mechanism like dispersion, absorption or reflections are not considered in this study.
If we consider that some average reaction temperature is achieved in both systems, then we could conclude that along the reactor an homogeneous reaction temperature distribution is obtained when it operates closely to steady state conditions, being this the best expected scenario. However, when the fluidized reactor is just starting to operate the temperature distribution of the medium is not homogeneous everywhere, affecting evidently the distribution of reacted particles. In fact at these operating conditions, the particles close to the walls of the first analyzed reactor will react faster than the ones placed at the center; and for the second type, the particles placed close to the reactor entrance or exit, will react more rapidly.

Diffusion-Convection Model
This problem is governed by Equation 7 without the source term Φ, leading to Equation A.1: x z t x z t Pe t dZ For the diffusion problem the convective term χPe∂θ/∂z is just dropped in the last equation. In order to find analytical solutions to both problems the following set of boundary conditions are established Equation A.2 and A.3: And the initial condition Equation A.4: Were-write last equations in terms of the new variable Θ(x, z, t) = θ(x, z, t)-1, the resulting equations are Equation A.5: With boundary condition Equation A.6 and A.7: And initial condition Equation A.8: The last problem can be solved by using the variables separation technique. In this sense the new temperature variable Θ(x, z, t), is expressed by Equation A.9: With boundary conditions Equation A.10 and A.11: Replacing Θ(x, z, t), in Equation A.15 and rearranging the terms we obtain Equation A.12: The last equation is only valid if both sides of equality are equal to a constant, for convenience, named as. The resulting set of Equations is A.13: And Equations (A.14): Similarly, the last equation is only valid if both sides of the equality are equal to a constant, also for convenience, named asµ. The resulting set of Equation A.15: where, λ and µ will be determined lines below.
where, n is an integer number and λ = n 2 π 2 /4. In a similar way, after applying the conditions in Equation where, the coefficient Equation A.23:

Diffusion-Convection Model with Heat Source
x z t x z t Pe dZ x z t x z t x z dX dZ For the case of pure diffusion, the convective term χPe∂θ/∂z is just dropped in the last equation. In order to find analytical solutions to both problems we assumed that the heat source term Φ is just a function of x. We also assumed the following set of boundary conditions Equation