- Email: [email protected]
On mathematical modeling of solar hydrogen production in monolithic reactors
Computers and Chemical Engineering 35 (2011) 1915–1922
Contents lists available at ScienceDirect
Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng
On mathematical modeling of solar hydrogen production in monolithic reactors M. Kostoglou a , C.P. Lekkos b , A.G. Konstandopoulos b,c,∗ a
Department of Chemistry, Aristotle University of Thessaloniki, Greece Aerosol and Particle Technology Laboratory (APTL), Chemical Process Engineering Research Institute (CPERI), Centre for Research and Technology Hellas (CERTH), Greece c Department of Chemical Engineering, Aristotle University of Thessaloniki, Greece b
a r t i c l e
i n f o
Article history: Received 17 November 2010 Received in revised form 22 March 2011 Accepted 23 March 2011 Available online 30 March 2011 Keywords: Solar hydrogen Water splitting Mathematical modeling Monolithic reactors
a b s t r a c t The hydrogen production based on the exploitation of the water splitting reaction and of solar energy on monolith type reactors is a novel and very promising process. Several trade offs between reaction efficiency and material stability appear during the operation of the process rendering necessary advanced optimization and control using appropriate mathematical models. The complexity of the geometric and operating features of the process prevents from the use of a general model for the optimization studies so a model splitting in several size scales procedure is proposed. It is shown that in case of silicon carbide as monolith material the one dimensional single channel model is the most appropriate to be used for the process operation optimization studies. The above general principles are described in detail in the present paper to shed light to the particular modeling problem with the emphasis given in the work performed in our laboratory. © 2011 Published by Elsevier Ltd.
1. Introduction Great efforts are concentrated worldwide on the conversion of solar energy into “solar fuels” with hydrogen being among the best candidates and water being the primary candidate raw material. Since direct one-step thermal dissociation of water (thermolysis) is a reaction requiring extremely high temperatures (>2200 ◦ C) for obtaining a significant dissociation degree, current research has shifted towards the so-called thermochemical water splitting cycles. These are “multi-step” processes where water is decomposed into hydrogen and oxygen via two or more chemical reactions taking place at much lower temperatures, that form a closed cycle wherein the overall reaction is: H2 O → H2 + ½ O2 and all other intermediate reactants and products are recycled (Kodama & Gokon, 2007; Steinfeld, Sanders, & Palumbo, 1999). The monolithic honeycomb solar reactor has been proposed for performing redox pair-based cycles for the production of hydrogen from the splitting of water using solar energy (Agrafiotis et al., 2005; Roeb et al., 2006) inspired by the well-known automobile catalytic converters. As steam passes through the solar reactor, the coating material (usually the lower-valence oxide of a metal exhibiting
∗ Corresponding author at: Aerosol and Particle Technology Laboratory (APTL), Chemical Process Engineering Research Institute (CPERI), Centre for Research and Technology Hellas (CERTH), 6th km Charilaou – Thermi Road, Thermi 57001, Thessaloniki, Greece. Tel.: +30 2310498192. E-mail addresses: [email protected] (M. Kostoglou), [email protected] (C.P. Lekkos), [email protected] (A.G. Konstandopoulos). 0098-1354/$ – see front matter © 2011 Published by Elsevier Ltd. doi:10.1016/j.compchemeng.2011.03.019
multiple oxidation states) splits water vapor by “trapping” its oxygen and leaving in the effluent gas stream pure hydrogen. In a subsequent step the oxygen-“trapping” coating is regenerated by increasing the amount of solar heat absorbed by the reactor and hence a cyclic operation is established. The most challenging step of the process is the regeneration step which is feasible in high temperature. For lower temperatures the regeneration reaction is not thermodynamically promoted whereas for higher temperature liquid phase appears, leading to the loss of activity of redox material and eventually to its inactivation (Allendorf, Diver, Siegel, & Miller, 2008). So the operation of the process requires a careful temperature manipulation with alternating cycles of “low temperature” (up to 1100 ◦ C) and “high temperature” (up to 1500 ◦ C). Depending on the particular redox material used this “high temperature” must be in a relatively narrow range in order to the process being efficient. Lower temperatures eliminate regeneration and higher temperatures completely destroy the redox material. The temperature handling is achieved through opening and closing the solar heat flux at the front of the monolith. The main problem is that the monolith has its thermal capacity and its internal heat transfer limitations rendering temperature control a difficult task. The temperature field in the monolith is non-uniform and transient and the only control variable for the temperature is the degree of solar radiation absorption. Another important control variable of the process is the composition of the inlet gas stream. Although the purpose is to use pure steam as working gas during the water splitting step, in practice, an inert gas is used restricting the fraction of steam to less than 20%. The inert gas is also used at the regeneration step to purge the
1916
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
produced oxygen. It is important to keep the two reactions (water splitting and regeneration) completely separated in time because the separation of hydrogen from oxygen is a difficult task. The inert gas can be nitrogen but preferably argon or helium to exclude the possibility of NOx formation during the regeneration step. So the inlet stream composition has two states. The control variable here is the time of switching from the one state to the other in respect to the solar radiation absorption cycle. The proper operation of the process needs a careful manipulation of the two control variables and in order to choose the appropriate operation scenario, a mathematical model of the process is necessary. The model needs to be accurate enough to predict local temperatures since they are crucial to ensure the activity of the redox material. A very specific optimization problem is set up for the particular problem at hand; the variable to be optimized is the average production rate of hydrogen which depends not only on the efficiency of the water splitting stage but on the ratio of the duration of the two operation stages. The simulation of the operation of monolithic reactors is in general a very demanding task due to the multiplicity of size and time scales to be accounted for. In particular regarding the size scales, the wall size scale, the channel size scale and the complete monolith size scale must be considered. Regarding time scales, the reaction time scale, the gas residence time scale and the externally imposed periodic operation time scale have to be resolved by the computational model. Of course, the problem of modeling transient operation of monolith is not new and it has been studied extensively in the literature, mainly in the context of catalytic automotive converters (Dubien, Schweich, Mabilon, Martin, & Prigent, 1998; Siemund, Leclerc, Schweich, Prigent, & Castagna, 1996). A recent state of the art review on available techniques for the simulation of automotive converters can be found in Mladenov, Koop, Tischer, and Deutschmann (2010). Nevertheless, the particular problem at hand has some distinct features that make it different than the automotive converters one. At first the much higher temperatures needed for the two reactions bring into play the radiation heat transfer mechanism which can be ignored for temperatures less than 600 ◦ C met in automotive converters (Young & Finlayson, 1976). In addition the large temperature differences developed between wall and gas in the channel set to question the perimeter averaged equations for the channels and the use of heat transfer coefficients. The issue of modeling of monolith reactors at large temperatures has also been considered for the problem of catalytic combustion (e.g. methane combustion; Tischer & Deutschmann, 2005; Tischer, Correa, & Deutschmann, 2001). The Fick law used for automotive converters considering the low concentration of reacting gases must be replaced by the Stefan Maxwell equations (Deen, 1998). Finally, the periodic operation introduces an additional time scale and the requirement of large time simulations. Given that the whole idea of producing hydrogen from steam in monolithic converters heated by the sun is new, only few modeling approaches can be found in the literature. A single perimeter averaged model can be found in Agrafiotis, Pagkoura, Lorenzou, Kostoglou, and Konstandopoulos (2007). A two dimensional pseudo-homogeneous model of the whole monolith using the simulation library Modelica can be found in Roeb et al. (2008) and Dersch, Mathijssen, Roeb, and Sattler (2006). In the present work a discussion of the modeling aspects of the recently proposed design for hydrogen production, will be presented with emphasis on the several levels of models needed to be developed in order to cover specific needs. The silicon carbide as monolith material shows distinct advantages (regarding its thermal properties) over traditional monolith materials as cordierite so it is the chosen material for the hydrogen production monoliths. The total efficiency of the process is assessed by two numbers. The first is the traditional degree of conversion of the vapour to hydrogen during the water
splitting step (reaction efficiency) and the second the fraction of cycle time spent in water splitting step (cycle efficiency). The structure of the present work is the following: At first the modeling problem is described following a classification per size scale. The corresponding governing equations are described in Appendix A to facilitate the reading of the main text. Then, indicative results are shown and discussed in order to confirm some of the statements of the previous section and to stress several aspects of the process. 2. Model development There are several aspects that must be included in a model of the solar hydrogen production process. In general, there are three coexisting size scales which must be taken into account. The first scale (level) refers to phenomena occurring on the monolith walls, the second level to phenomena occurring along the monolith channels and the third level to phenomena occurring at the macroscopic scale of the monolith. Submodels for each level must be developed and properly combined in order to derive a model for the entire process. In the following the three levels are discussed one by one. 2.1. Level I, monolith wall scale The study of the particular subject is still at first stages so advanced issues as the structural details of the wall/redox material and their influence on the reaction rates (further to some simple extensions as the shrinking core model) have not discussed yet. It is noted that structural details have been considered for similar problems met in case of monolith converters (Koci, Kubicek, & Marek, 2004) or diesel particulate filters (Kostoglou & Konstandopoulos, 2005). So ignoring the structure of the redox material, the model of this scale degenerates to expressions for the kinetics of the reactions. In a very simplified presentation the occurred reactions can be described as Water splitting (exothermic reaction) M + H2 O → MO + H2
(1a)
Regeneration (endothermic reaction) MO → M +
1 O2 2
(1b)
where M and MO are the reduced and oxidized forms of the redox material respectively. The oxidized form MO belongs to the family of ferrites, i.e. it has the chemical structure XFe2 O4 where X can be Fe, Mn, Co, Ni, Zn. Extensive thermodynamic calculations reveal the potential of the ferrites to work properly at the above cycle (Roeb, Gathmann, Neises, Sattler, & Pitz-Paal, 2009). The crucial step is the regeneration (or as alternatively called thermal reduction) reaction. The product of the regeneration reaction should be one that allows the continuation of the thermochemical cycles and it should be taken in temperatures lower than those that can cause destruction of the material. For example, the regeneration produces the required material for X = Fe, Zn so high temperatures for which the material suffers sintering (X = Fe) or evaporates (Zn). The regeneration reaction gives the required product only in a narrow window of temperature. For higher temperatures, liquid instead of solid product appears reducing the efficiency of the process and leading to its termination after a small number of cycles. Further to the thermodynamics, the kinetics of the reactions (1) is very important. The water splitting reaction may be actually consisted of several elementary steps. What is needed for the model development is not the exact reaction mechanism but the simplest
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
overall reaction rate expression that can describe the experimental data. An empirical water splitting kinetics scheme is proposed by Tsuji, Togawa, Wada, Sano, & Tamaura (1995) according to which the reaction rate depends only on the partial pressures of vapor and hydrogen and not on the state of the redox material. A much simpler rate expression containing a linear dependence on the vapor concentration is proposed by Chen and Hogan (2009). A fairly general reaction rate scheme has the form (Rs is the water splitting reaction rate in moles of produced hydrogen per unit area of redox material per second): Rs =
ks (1 − y)C 1 + kH C
(2)
where ks is a temperature dependent rate constant, is the maximum oxygen storage capacity of the redox material (total moles H2 produced per redox material surface area unit) and it also depends on the temperature and kH is a temperature dependent constant related to surface adsorption–desorption process. Finally, C is the vapor concentration and y is the instantaneous actual storage normalized by the maximum one and it takes values between 0 and 1. The quantity (1 − y) is proportional to the number of empty sites available for oxygen storage on the material. This rate expression is the outcome of combined theoretical considerations and experimental observations. Regarding regeneration, a constant rate has been proposed (Chen & Hogan, 2009) but the most reasonable choice is the first order dependence on the adsorbed oxygen i.e. Rr = kr y
(3)
where the reaction rate constant kr exhibits an Arrhenius temperature dependence. The water splitting and regeneration reaction constants are computed from experimental data taken in packed bed reactors with the redox material in the form of powders (Zygogianni et al., 2009). 2.2. Level II, monolith channel scale The transient operation of the monolith has in general a time scale large compared to the residence time of the gas in the monolith so a quasi-steady approximation for the flow of the gas can be considered. The one dimensional (1D) model of this scale includes the quasi-steady momentum, heat and mass balances for the gas stream (see Appendix A.1 for details). The diffusion of water vapor towards the wall has been modeled in Agrafiotis, Lorentzou, Pagkoura, Kostoglou, and Konstandopoulos (2006) using the Fick law that allows an analytical solution to the vapor mass balance in the gas phase. The use of Fick law is justified in those cases since the focus of the work was on modeling experimental campaigns with vapor molar fraction less than 10% in the feed stream. In general, considering larger fraction of vapor, the diffusion process must be described by Stefan Maxwell equations. A need arise here to develop perimeter averaged 1-dimensional models based on Stefan Maxwell equations instead on the Fick law. The key scope is to include the generality of Stefan Maxwell equations keeping the computational cost at minimum. Another pitfall of the existing 1D model is the ignorance of the radiation heat transfer mechanism. The physical properties of the gas are computed at the local cross sectional average value of the temperature but although it is a standard practice for the existing monolith applications it remains to be confirmed for the high temperature range met in the present application. Even with the above pitfalls the 1D model was able to predict the trends of the experimental results for hydrogen production cycles. An important computational aspect revealed by the 1D model is the stiffness of the set of equations regarding their integration in the time domain. Whereas in case of automotive converters the stiffness is due to the extremely high temperature dependence
1917
of the Langmuir–Hinshelwood expressions for the kinetic of the oxidation reactions (Oh & Cavendish, 1982), in the present case the stiffness originates from the heat transfer problem itself. The solar radiation enters the monolith as a flux boundary condition on the solid face of the monolith and the heat is moving slowly into the solid wall as a steep front. Unfortunately, this stiffness is inherent to the process and it is an obstacle in the use of higher level models to simulate the operation of the hydrogen operation unit. In order to examine the validity of the 1D model and to adjust it to accommodate the additional complexities of the actual problem (Stefan Maxwell mass transfer, radiation intra-channel heat transfer and high cross-sectional temperature non-uniformities), a computational fluid dynamics (CFD) oriented single channel model is developed (see Appendix A.2 for its description). The channel interior and the channel walls are discretized and the continuity equation, the momentum balances and the heat and mass conservation equations employing Stefan Maxwell diffusive flux relations, are solved in the gas phase combined with the heat conservation equation, the reaction kinetic rates and the equation for the accumulation of the oxygen in the solid phase. The computational cost of the three dimensional (3D) single channel model is very high so the effects of the 3rd dimension and of the radiation are studied carefully during a heat up of the monolith with solar radiation, up to the establishment of the steady state. Using the two dimensional (2D) equivalent problem with properly transformed variables it can be shown that the information expected from the 3D model can be also extracted using the 2D model. The CFD code admits detailed inclusion of the surface to surface radiation model. The analysis of the heat up showed that radiation acts as an additional thermal conductivity of the solid wall. This conductivity is highly depended on the temperature (monotonically increasing function of temperature). For materials with low thermal conductivity as the cordierite, the effect of the radiation on heat transfer in the monolith is important and it cannot be ignored. On the contrary the high thermal conductivity of the silicon carbide makes the contribution of the radiation negligible up to 1500 ◦ C so it can be ignored for the simulation of the cycling operation of the monolithic reactor reducing decisively the computational requirements. 2.3. Level III, entire reactor scale This is the macroscopic scale of the monolith. The need for models of this scale is due to the spatial non-uniformities of the inlet conditions or to boundary heat loses. For example, in the case of the automotive converters and diesel particulate filters (Kostoglou, Housiada, & Konstandopoulos, 2003) the non-uniformities are met mainly at the inlet flow distribution, at the inlet gas temperature distribution and they are created by the heat loses from the monolith to the environment. The above non-uniformity sources are of minor significance for the hydrogen production monolith reactor for which the primary cause of non-uniformity is the radiation strength distribution. Although, a perfectly flat heliostat would produce an image on the receiver of the size of the heliostat, for most applications each mirror segment on a heliostat is concaved slightly and canted towards the receiver. Mirror surface waviness, the gross curvature error of each heliostat, the errors associated with accurate positioning of each mirror segment and the effects of differential thermal growth tend to increase the image error. All of these add up optically to produce a higher flux density at the aim point which has a distribution pattern of a Gaussian shape typically used in simulation studies (Stine & Geyer, 2001). The detailed way to simulate the entire monolith is to discretize all the channels one by one but this approach is computationally out of question. A more reasonable approach is to simulate all channels using the 1D perimeter averaged model and allow communication between the channels through the heat transferred among their
1918
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
Fig. 1. Effect of radiation on the wall temperature profile evolution during monolith heating.
walls. Considering that a monolith typically contains a few thousands of channels this model is also computationally intractable. The usual approach for monolith scale simulations is to exploit the small value of the ratio between the channel side size and the monolith face area size and to homogenize the discrete channel set of equations to get a continuum model. The resulting homogeneous problem in which the channel level structure has been lost but it is indirectly conveyed through the appropriate transport coefficients can be solved by CFD codes as porous media. Many commercial CFD codes assume porous media to be always in thermal equilibrium between gas and solid phase. This problem is overcome by the so called dual grid approach according to which a solid phase and a porous medium are assumed to occupy the same space and connected to each other through the appropriate transport relations. A further discussion on the models of this level can be found in Appendix A.3. 3. Indicative results and discussion There are two modes of operation of the monolith: the first is the relative long heating period required for the cold monolith to reach the reaction temperature. During this period only inert gas is fed to the monolith and the only occurring phenomenon is heat transfer. The second mode is the periodic operation during which the reactions occur and the hydrogen is produced. The 3D single channel problem with radiation is too expensive computationally to be used for the simulation of the hydrogen production cycles. Comparison between the simple heating (first mode) results taken for 2D (parallel plates) and 3D (square cross section) channels revealed only quantitative differences that can be accommodated by modifying appropriately the heat transfer coefficient. So considering that the purpose of the CFD model is the adjustment of the 1D perimeter averaged single channel model, the 2D channel geometry can be employed. The inclusion of the radiation also requires significant computational effort. The wall surface temperature profile along the channel for typical heating conditions is shown in Fig. 1 in the absence and presence of radiation. The high thermal conductivity of the silicon carbide leads to a relative uniform temperature profile which does not give space to radiation for further smoothing of this profile. The effect of the radiation is negligible so radiation can be completely ignored for monoliths made by silicon carbide. The production cycles are simulated using the 2-D geometry without radiation. The simulated scenario is the following: In order to converge fast to a limit cycle, a steady state operation at 1660 K with 100% nitrogen as gas phase and all the redox material in the reduced state is assumed as initial condition. Then the solar heat flux decreases and when the temperature goes down from the 1473 K the feed changes to 80% nitrogen and 20% water vapor.
Fig. 2. Channel wall average temperature and gas outlet temperature during the periodic operation of the monolith.
When the temperature reaches the 1073 K the solar heat flux goes again to its initial value and at 1473 K the feed change again to 100% nitrogen. The convergence to a periodic steady state occurs at the second cycle. The evolution of the average wall temperature and of the gas temperature at the outlet of the monolith is presented in Fig. 2. The second cycle is the limiting one and the system has come to a periodic steady state. There is a delay between the gas and the average solid temperature. The outlet gas is in thermal equilibrium with the outflow edge of the monolith wall which always exhibits a delay with respect to the average wall temperature since the heat source is located at the inlet edge. The evolution of the concentration of the important species in the outflow is shown in Fig. 3. During the water splitting stage the vapor is completely transformed to hydrogen. When most of the redox material is in the oxidized state the rate of the reaction is reduced and a part of the vapor remains unreacted. The reaction efficiency of the water splitting process is 83% and the cycle efficiency (duration of water splitting step to total cycle duration) is 54% in the present case. The evolution of the fraction of the redox material in the reduced and oxidized states is shown in Fig. 4. The present case corresponds to an efficient operation scenario since the water splitting stage stops at oxidized state fraction y close to 0.8 and starts at y close to 0.05. The surface wall temperature profiles along the channel at several moments during the cycle are shown in Fig. 5. As it is obvious, the temperature profiles along the channel are highly non-uniform due to the combination of the strongly localized heating (entrance of the channel) and the need for periodic operation. The nonuniform conditions along the filter combined to the very strict requirements on the local temperature range (especially during the regeneration step) renders the selection and optimization of the proper operation strategy, a very demanding task. The evolution of
Fig. 3. Chemical species concentration at the outlet gas stream versus time.
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
Fig. 4. Evolution of the (complementary) fractions of oxidized (MO) and reduced (M) redox material.
1919
Fig. 7. Contour plot of the hydrogen mass fraction in the channel of the monolith at t = 720 s.
Fig. 8. Temperature (K) contour plot in the channel of the monolith at t = 720 s.
Fig. 5. Typical temperature profiles along the channel wall indicating the extent of temperature non-uniformity.
the fraction y of the oxidized redox material along the channel is shown in Fig. 6. Both oxidation and regeneration proceed from the inlet to the outlet creating highly non-uniform profiles of y. The distribution of the hydrogen in the channel during the water splitting stage is shown in Fig. 7 (at t = 720 s). The occurrence of the reaction on the wall leads to large hydrogen concentrations close to the wall. The vapor is completely consumed along the flow and the already produced hydrogen is transferred from the wall to the whole channel cross section by diffusion. The species concentration is reduced along the flow due to the gas thermal expansion as the temperature increases. The temperature contours in the channel at the same time (t = 720 s), are shown in Fig. 8. The gas is losing fast its thermal content due to its low thermal capacity and it is adjusted
Fig. 6. Evolution of profiles of the fraction y of oxidized redox material along the channel.
to follow the wall temperature. The wall temperature is uniform across the wall and varies only along the channel. The gas temperature follows the longitudal variation of the wall temperature profile retaining a self-similar transverse profile that corresponds to the asymptotic Nusselt number. The whole structure of the temperature field in Fig. 8 indicates that a 1D approach can accurately predict the heat transfer dynamics in the channel. The full scale monolith model is needed in order to study the effect of heat source non-uniformities and macroscopic geometry. A typical industrial scale setup of the monolith was studied during its heating by a heat source of Gaussian type with its maximum at the center line of the monolith and a close to zero value of flux at the circumference. Several steady state results indicating the effect of source non-uniformity on the temperature non-uniformities are presented. The temperature contours and the gas phase streamlines are shown in Fig. 9. At this temperature scale transverse nonuniformities of the temperature are not visible. The small thermal capacity of the gas in the mixing chamber result in its preheating before entering the monolith so the gas stream enters into the monolith in thermal equilibrium. The effect of the macroscopic equipment structure must be considered by channel level models that typically assume different temperature between monolith and inlet gas. The temperature profiles on three types of commercial monoliths (with 90, 180 and 270 channels per square inches – cpsi) are shown in Fig. 10. The large heat source non-uniformity has been drastically smeared-out by the high conductivity of the material leading to small temperature non-uniformities. As the size of the
Fig. 9. Temperature (K) contours (a) and gas flow streamlines (b) for a heated by non-uniform solar heat source monolith at steady state.
1920
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
Fig. 10. Steady state temperature (K) contours for monoliths (with different channel densities) heated by non-uniform solar heat source.
Fig. 12. Schematic and variables of single channel.
Appendix A. A.1. 1D single channel model
Fig. 11. Wall temperature (K) for channels at the center and at the circumference of monoliths (with different channel densities) heated by non-uniform solar heat source.
channel decreases the effective thermal conductivity of the monolith increases leading to smaller temperature non-uniformities as shown in Fig. 10. The temperature profiles along the channel for channels being at the center and at the circumference of the monolith are presented in Fig. 11. The non-uniformities at the face cross section are eliminated by the conductance along the channel as heat “diffuses” from the hotter central area to the colder circumference of the monolith. Conclusively, one could say that at least for silicon carbide the temperature non-uniformities created by the solar flux non-uniformities are small enough to be ignored in operation strategy optimization studies.
The equations that describe the single channel dynamics are presented in the following (where ˛ is the side size of the square channel, L the length of monolith and w the thickness of the wall of the monolith – see Fig. 12). The unknowns of the model are the evolved spatial (along the channel) profiles of gas phase composition and temperature, of wall temperature and of state y of the redox material. Mass conservation equations are preferably written in terms of molar flux of species instead of concentrations in order to avoid complications arising by the gas density variation (A is the redox material surface area per filter wall surface area):
A.1.1. Gas phase Species conservation: (i) H2 O (vapor) df1 = −4˛kD (c1 − c1w ) dz
(A1)
The flux to the wall must be equal to the reaction rate: kD (cH − cHw ) = ARs
(A2)
(ii) inert gas 4. Conclusions The need to optimize the operation of hydrogen production from heated monoliths process asks for advanced modeling of the process. Unfortunately, the disparate size and time scales appeared in the problem do not admit the use of a general purpose simulation tool for the demanding optimization study of the process. Different models must be used for different purposes. In this view, the three-dimensional complete monolith and housing model is necessary to study (under simple operations conditions as pure heating or steady state) the effect of heat flux or flow non-uniformity to the uniformity of the conditions at the outlet of the monolith. The detailed CFD single channel model is needed to properly adjust and validate 1D perimeter averaged single channel model under typical operation conditions. Finally, the properly adjusted 1D channel model is the appropriate scale model to be used for the extensive computations needed in view of optimization of hydrogen production unit operation. Regarding monolith material, it was shown that the high thermal conductivity of the silicon carbide simplifies considerably the modeling requirements by eliminating transverse temperature gradients and allowing ignorance of the intrachannel radiation heat transfer. The computational requirements are much higher for low conductivity materials like cordierite in which steep temperature fronts appear making necessary the consideration of transverse non-uniformity and intrachannel radiation by the model.
df2 =0 dz
(A3)
(iii) H2 df3 = 4˛ARs dz
(A4)
(iv) O2 df4 = 4˛ARr dz
(A5)
Energy conservation:
4
fi cpi
dT = 4˛kh (Tw − T ) dz
(A6)
i=1
A.1.2. Solid phase Surface coverage evolution:
dy = ARs − ARr dt
(A7)
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
Energy conservation: 2s cpw w˛
∂Tw = 2w w˛ ∂t
∂2 T ∂z 2
+ 4˛kh (T − Tw )
+ 4˛A(Hs Rs + Hr Rr )
(A8)
The concentrations in the gas phase is computed as P ci = 4 R f gT i=1 i
(A9)
In the above equations cpi is the molar specific heat capacity of the gaseous species i, cpw the specific heat capacity of the channel wall, , s the thermal conductivity and density of the wall, H are the reaction enthalpies, kh the channel to wall heat transfer coefficient and kD the channel to wall mass transfer coefficient. The molar fluxes of vapor, inert gas, hydrogen and oxygen are denoted as f1 , f2 , f3 and f4 respectively. Finally P, T denote the local pressure and temperature and Rg the gas constant. The conductivity of the porous wall is computed using typical expressions for the conductivity of composite media and the temperature dependent intrinsic conductivity of the silicon carbide. In case of siliconized silicon carbide being the monolith material, the solid phase intrinsic conductivity is computed through a weighted average of silica and silicon carbide. The mass transfer coefficient kD can be modified to account for the effect of Stefan Maxwell equations without increasing the mathematical complexity of the problem. Instead of the ternary mixture of water vapour, hydrogen and inert gas a binary system of vapour diffused in a mixture of hydrogen and inert gas is considered. Extensions of Fick’s law for binary mixtures and for diffusion in mixture of gases can be combined to give a concentration dependent diffusion coefficient which through the Sherwood number relation for square channel (asymptotic Sherwood number equal to 3) leads to the coefficient kD . The resulting relation for the effective diffusivity of vapor in the mixture of inert gas and hydrogen is Def = (xH /DwH + xin /Dwin )−1 where xH , xin are the local molar fractions of hydrogen and inert gas respectively, and DwH , Dwin are the binary diffusivities of vapor in hydrogen and inert gas respectively. A.1.3. Boundary conditions The boundary conditions for the monolith are as follows: fi (0, t) = fio (t)
(A10a)
T (0, t) = To (t)
(A10b) at z = 0
∂Tw = 0 at z = L ∂z
(A10c) (A10d)
where Q(t) is the solar power adsorbed by the monolith per channel. A.1.4. Initial condition Tw (x, 0) = Two
coupling between gas and solid phase is achieved through continuity of temperature, heat flux and reacting species fluxes. The surface to surface radiation model is considered since the gas filling the channels does not absorb radiation. Finally the temperature dependence of the gas properties is taken into account. A.3. 3D model
fi
∂Tw Q (t) − = ∂z (2w˛ + w2 )
1921
(A10e)
A.2. CFD-based detailed single channel model Only a description of the mathematical problem for this level model will be given here to avoid copying of equations described in the manual of a commercial CFD code (specifically the code FLUENT has been used here). The governing equations are the continuity equation, the momentum conservation equations, the Stefan Maxwell relations for the diffusive fluxes, the gas phase energy conservation equation and the solid phase thermal conservation equation. The surface reaction feature of the code is used for the consideration of water splitting and regeneration reactions. The
The three dimensional model results from a homogenization procedure during which the geometrical details of the monolith are lost and the solid and gas phases are considered as interpenetrating continua. The governing equations are the same with those of the 1dimensional single channel model but with two differences: (i) Two additional spatial independent variables are needed (x, y) rendering all the depended variables functions of x, y, z, t. (ii) A mesoscopic conduction term has to been added. The effective conductivity tensor has two components one for the axial direction and one for the transverse direction. Derivation of these components for monoliths with square channels can be found elsewhere (Groppi & Tronconi, 1996; Konstandopoulos, Kostoglou, Vlachos, & Kladopoulou, 2007). The resulting system of equations can be solved by assuming representative channels on the monolith cross-section, solving them as 1D problems, find the heat sources for each channel and update the solid phase heat transfer equation with these sources. Alternatively, the CFD code FLUENT can be used through the so called dual grid approach considering a porous medium and a compact solid properly connected to each other. The reactor scale results presented here were taken employed this approach. References Agrafiotis, C., Roeb, M., Konstandopoulos, A. G., Nalbandian, L., Zaspalis, V. T., & Sattler, C. (2005). Solar water splitting for hydrogen production with monolithic reactors. Journal of Solar Energy Engineering, 79, 409–421. Agrafiotis, C., Lorentzou, S., Pagkoura, C., Kostoglou, M., & Konstandopoulos, A. G. (2006). Advanced monolithic reactors for hydrogen generation from solar water splitting. In Seville: Proceedings of solarPACES 13th international symposium on concentrated solar power and chemical energy technologies 2006 Seville, Spain, June 20–23. Agrafiotis, C., Pagkoura, C., Lorenzou, C., Kostoglou, M., & Konstandopoulos, A. G. (2007). Hydrogen production in solar reactors. Catalysis Today, 127, 265–277. Allendorf, M. D., Diver, R. B., Siegel, N. P., & Miller, J. E. (2008). Two step water splitting using mixed-metal ferrites: Thermodynamic analysis and characterization of synthesized materials. Energy & Fuels, 22, 4115–4124. Chen, K. S., & Hogan, R. E. (2009). A two phase model for solar thermochemical water splitting with FeO/Fe3 O4 . In Proceedings of ES2009. 2009 ASME 3rd international conference of energy sustainability San Francisco, California, July 19–23. Deen, W. M. (1998). Analysis of transport phenomena. Oxford University Press: New York. Dersch, J., Mathijssen, A., Roeb, M., & Sattler, C. (2006). Modelling of a solar thermal reactor for hydrogen generation. Modelica. Dubien, C., Schweich, D., Mabilon, G., Martin, B., & Prigent, M. (1998). Three-way catalytic converter modelling: Fast- and slow-oxidizing hydrocarbons, inhibiting species, and steam-reforming reaction. Chemical Engineering Science, 53, 471–481. Groppi, G., & Tronconi, E. (1996). Continuous vs. discrete models of non-adiabatic monolith catalysts. AIChE Journal, 42, 2382–2387. Koci, P., Kubicek, M., & Marek, M. (2004). Modeling of three-way-catalyst monolith converters with microkinetics and diffusion in the washcoat. Industrial & Engineering Chemistry Research, 43, 4503–4510. Kodama, T., & Gokon, N. (2007). Thermochemical cycles for high-temperature solar hydrogen production. Chemical Reviews, 107, 4048–4077. Konstandopoulos, A. G., Kostoglou, M., Vlachos, N., & Kladopoulou, E. (2007). Advances in the science and technology of diesel particulate filter simulation. In G. B. Marin (Ed.), Advances in chemical engineering, Special issue on automotive emission control (pp. 213–275). Academic Press. Kostoglou, M., & Konstandopoulos, A. G. (2005). Effect of soot layer microstructure on diesel particulate filter regeneration. AIChE Journal, 51, 2534–2546. Kostoglou, M., Housiada, P., & Konstandopoulos, A. G. (2003). Multi-channel simulation of regeneration in honeycomb monolithic diesel particulate filters. Chemical Engineering Science, 58, 3273–3283. Mladenov, N., Koop, J., Tischer, S., & Deutschmann, O. (2010). Modeling of transport and chemistry in channel flows of automotive catalytic converters. Chemical Engineering Science, 65, 812–826. Oh, S. H., & Cavendish, J. C. (1982). Transients of monolithic catalytic converters: Response to step changes in feedstream temperature as related to controlling
1922
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
automobile emissions. Industrial & Engineering Chemistry Product Research and Development, 21, 29–37. Roeb, M., Sattler, C., Klüser, R., Monnerie, N., de Oliveira, L., & Konstandopoulos, A. G. (2006). Solar hydrogen production by a two-step cycle based on mixed iron oxides. Journal of Solar Energy Engineering, 128, 125–133. Roeb, M., Neises, M., Säck, J. P., Rietbrock, P., Monnerie, N., & Dersch, J. (2008). Operational strategy of a two-step thermochemical process for solar hydrogen production. International Journal of Hydrogen Energy, doi:10.1016/j.ijhydene.2008.08.049 Roeb, M., Gathmann, N., Neises, M., Sattler, C., & Pitz-Paal, R. (2009). Thermodynamic analysis of the two-step solar water splitting with mixes iron oxides. International Journal of Energy Research, doi:10.1002/er.1513 Siemund, S., Leclerc, J. P., Schweich, D., Prigent, M., & Castagna, F. (1996). Three-way monolithic converter: Simulations versus experiments. Chemical Engineering Science, 51, 3709–3720. Steinfeld, A., Sanders, S., & Palumbo, R. (1999). Design aspects of solar thermochemical engineering – A case study: Two-step water-splitting cycle using the Fe3 O4 /FeO redox system. Solar Energy, 65, 43–53.
Stine, W. B., & Geyer, M. (2001). Power from the sun (electronic edition). Tischer, S., & Deutschmann, O. (2005). Recent advances in numerical modeling of catalytic monolith reactors. Catalysis Today, 105, 407– 413. Tischer, S., Correa, C., & Deutschmann, O. (2001). Transient three-dimensional simulations of a catalytic combustion monolith using detailed models for heterogeneous and homogeneous reactions and transport phenomena. Catalysis Today, 69, 57–62. Tsuji, M., Togawa, T., Wada, Y., Sano, T., & Tamaura, Y. (1995). Kinetic study of the formation of cation-excess magnetite. Journal of the Chemical Society – Faraday Transactions, 91, 1533–1538. Young, L. C., & Finlayson, B. A. (1976). Mathematical-models of monolith catalytic-converter. Application to automobile exhaust. AIChE Journal, 22, 343– 353. Zygogianni, A., Lorentzou, S., Pagkoura, C., Lekkos, C., Kostoglou, M., & Agrafiotis, C. (2009). On the modeling of solar monolith reactors for hydrogen production. In Proceedings of SolarPACES 2009 conference Berlin.
Contents lists available at ScienceDirect
Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng
On mathematical modeling of solar hydrogen production in monolithic reactors M. Kostoglou a , C.P. Lekkos b , A.G. Konstandopoulos b,c,∗ a
Department of Chemistry, Aristotle University of Thessaloniki, Greece Aerosol and Particle Technology Laboratory (APTL), Chemical Process Engineering Research Institute (CPERI), Centre for Research and Technology Hellas (CERTH), Greece c Department of Chemical Engineering, Aristotle University of Thessaloniki, Greece b
a r t i c l e
i n f o
Article history: Received 17 November 2010 Received in revised form 22 March 2011 Accepted 23 March 2011 Available online 30 March 2011 Keywords: Solar hydrogen Water splitting Mathematical modeling Monolithic reactors
a b s t r a c t The hydrogen production based on the exploitation of the water splitting reaction and of solar energy on monolith type reactors is a novel and very promising process. Several trade offs between reaction efficiency and material stability appear during the operation of the process rendering necessary advanced optimization and control using appropriate mathematical models. The complexity of the geometric and operating features of the process prevents from the use of a general model for the optimization studies so a model splitting in several size scales procedure is proposed. It is shown that in case of silicon carbide as monolith material the one dimensional single channel model is the most appropriate to be used for the process operation optimization studies. The above general principles are described in detail in the present paper to shed light to the particular modeling problem with the emphasis given in the work performed in our laboratory. © 2011 Published by Elsevier Ltd.
1. Introduction Great efforts are concentrated worldwide on the conversion of solar energy into “solar fuels” with hydrogen being among the best candidates and water being the primary candidate raw material. Since direct one-step thermal dissociation of water (thermolysis) is a reaction requiring extremely high temperatures (>2200 ◦ C) for obtaining a significant dissociation degree, current research has shifted towards the so-called thermochemical water splitting cycles. These are “multi-step” processes where water is decomposed into hydrogen and oxygen via two or more chemical reactions taking place at much lower temperatures, that form a closed cycle wherein the overall reaction is: H2 O → H2 + ½ O2 and all other intermediate reactants and products are recycled (Kodama & Gokon, 2007; Steinfeld, Sanders, & Palumbo, 1999). The monolithic honeycomb solar reactor has been proposed for performing redox pair-based cycles for the production of hydrogen from the splitting of water using solar energy (Agrafiotis et al., 2005; Roeb et al., 2006) inspired by the well-known automobile catalytic converters. As steam passes through the solar reactor, the coating material (usually the lower-valence oxide of a metal exhibiting
∗ Corresponding author at: Aerosol and Particle Technology Laboratory (APTL), Chemical Process Engineering Research Institute (CPERI), Centre for Research and Technology Hellas (CERTH), 6th km Charilaou – Thermi Road, Thermi 57001, Thessaloniki, Greece. Tel.: +30 2310498192. E-mail addresses: [email protected] (M. Kostoglou), [email protected] (C.P. Lekkos), [email protected] (A.G. Konstandopoulos). 0098-1354/$ – see front matter © 2011 Published by Elsevier Ltd. doi:10.1016/j.compchemeng.2011.03.019
multiple oxidation states) splits water vapor by “trapping” its oxygen and leaving in the effluent gas stream pure hydrogen. In a subsequent step the oxygen-“trapping” coating is regenerated by increasing the amount of solar heat absorbed by the reactor and hence a cyclic operation is established. The most challenging step of the process is the regeneration step which is feasible in high temperature. For lower temperatures the regeneration reaction is not thermodynamically promoted whereas for higher temperature liquid phase appears, leading to the loss of activity of redox material and eventually to its inactivation (Allendorf, Diver, Siegel, & Miller, 2008). So the operation of the process requires a careful temperature manipulation with alternating cycles of “low temperature” (up to 1100 ◦ C) and “high temperature” (up to 1500 ◦ C). Depending on the particular redox material used this “high temperature” must be in a relatively narrow range in order to the process being efficient. Lower temperatures eliminate regeneration and higher temperatures completely destroy the redox material. The temperature handling is achieved through opening and closing the solar heat flux at the front of the monolith. The main problem is that the monolith has its thermal capacity and its internal heat transfer limitations rendering temperature control a difficult task. The temperature field in the monolith is non-uniform and transient and the only control variable for the temperature is the degree of solar radiation absorption. Another important control variable of the process is the composition of the inlet gas stream. Although the purpose is to use pure steam as working gas during the water splitting step, in practice, an inert gas is used restricting the fraction of steam to less than 20%. The inert gas is also used at the regeneration step to purge the
1916
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
produced oxygen. It is important to keep the two reactions (water splitting and regeneration) completely separated in time because the separation of hydrogen from oxygen is a difficult task. The inert gas can be nitrogen but preferably argon or helium to exclude the possibility of NOx formation during the regeneration step. So the inlet stream composition has two states. The control variable here is the time of switching from the one state to the other in respect to the solar radiation absorption cycle. The proper operation of the process needs a careful manipulation of the two control variables and in order to choose the appropriate operation scenario, a mathematical model of the process is necessary. The model needs to be accurate enough to predict local temperatures since they are crucial to ensure the activity of the redox material. A very specific optimization problem is set up for the particular problem at hand; the variable to be optimized is the average production rate of hydrogen which depends not only on the efficiency of the water splitting stage but on the ratio of the duration of the two operation stages. The simulation of the operation of monolithic reactors is in general a very demanding task due to the multiplicity of size and time scales to be accounted for. In particular regarding the size scales, the wall size scale, the channel size scale and the complete monolith size scale must be considered. Regarding time scales, the reaction time scale, the gas residence time scale and the externally imposed periodic operation time scale have to be resolved by the computational model. Of course, the problem of modeling transient operation of monolith is not new and it has been studied extensively in the literature, mainly in the context of catalytic automotive converters (Dubien, Schweich, Mabilon, Martin, & Prigent, 1998; Siemund, Leclerc, Schweich, Prigent, & Castagna, 1996). A recent state of the art review on available techniques for the simulation of automotive converters can be found in Mladenov, Koop, Tischer, and Deutschmann (2010). Nevertheless, the particular problem at hand has some distinct features that make it different than the automotive converters one. At first the much higher temperatures needed for the two reactions bring into play the radiation heat transfer mechanism which can be ignored for temperatures less than 600 ◦ C met in automotive converters (Young & Finlayson, 1976). In addition the large temperature differences developed between wall and gas in the channel set to question the perimeter averaged equations for the channels and the use of heat transfer coefficients. The issue of modeling of monolith reactors at large temperatures has also been considered for the problem of catalytic combustion (e.g. methane combustion; Tischer & Deutschmann, 2005; Tischer, Correa, & Deutschmann, 2001). The Fick law used for automotive converters considering the low concentration of reacting gases must be replaced by the Stefan Maxwell equations (Deen, 1998). Finally, the periodic operation introduces an additional time scale and the requirement of large time simulations. Given that the whole idea of producing hydrogen from steam in monolithic converters heated by the sun is new, only few modeling approaches can be found in the literature. A single perimeter averaged model can be found in Agrafiotis, Pagkoura, Lorenzou, Kostoglou, and Konstandopoulos (2007). A two dimensional pseudo-homogeneous model of the whole monolith using the simulation library Modelica can be found in Roeb et al. (2008) and Dersch, Mathijssen, Roeb, and Sattler (2006). In the present work a discussion of the modeling aspects of the recently proposed design for hydrogen production, will be presented with emphasis on the several levels of models needed to be developed in order to cover specific needs. The silicon carbide as monolith material shows distinct advantages (regarding its thermal properties) over traditional monolith materials as cordierite so it is the chosen material for the hydrogen production monoliths. The total efficiency of the process is assessed by two numbers. The first is the traditional degree of conversion of the vapour to hydrogen during the water
splitting step (reaction efficiency) and the second the fraction of cycle time spent in water splitting step (cycle efficiency). The structure of the present work is the following: At first the modeling problem is described following a classification per size scale. The corresponding governing equations are described in Appendix A to facilitate the reading of the main text. Then, indicative results are shown and discussed in order to confirm some of the statements of the previous section and to stress several aspects of the process. 2. Model development There are several aspects that must be included in a model of the solar hydrogen production process. In general, there are three coexisting size scales which must be taken into account. The first scale (level) refers to phenomena occurring on the monolith walls, the second level to phenomena occurring along the monolith channels and the third level to phenomena occurring at the macroscopic scale of the monolith. Submodels for each level must be developed and properly combined in order to derive a model for the entire process. In the following the three levels are discussed one by one. 2.1. Level I, monolith wall scale The study of the particular subject is still at first stages so advanced issues as the structural details of the wall/redox material and their influence on the reaction rates (further to some simple extensions as the shrinking core model) have not discussed yet. It is noted that structural details have been considered for similar problems met in case of monolith converters (Koci, Kubicek, & Marek, 2004) or diesel particulate filters (Kostoglou & Konstandopoulos, 2005). So ignoring the structure of the redox material, the model of this scale degenerates to expressions for the kinetics of the reactions. In a very simplified presentation the occurred reactions can be described as Water splitting (exothermic reaction) M + H2 O → MO + H2
(1a)
Regeneration (endothermic reaction) MO → M +
1 O2 2
(1b)
where M and MO are the reduced and oxidized forms of the redox material respectively. The oxidized form MO belongs to the family of ferrites, i.e. it has the chemical structure XFe2 O4 where X can be Fe, Mn, Co, Ni, Zn. Extensive thermodynamic calculations reveal the potential of the ferrites to work properly at the above cycle (Roeb, Gathmann, Neises, Sattler, & Pitz-Paal, 2009). The crucial step is the regeneration (or as alternatively called thermal reduction) reaction. The product of the regeneration reaction should be one that allows the continuation of the thermochemical cycles and it should be taken in temperatures lower than those that can cause destruction of the material. For example, the regeneration produces the required material for X = Fe, Zn so high temperatures for which the material suffers sintering (X = Fe) or evaporates (Zn). The regeneration reaction gives the required product only in a narrow window of temperature. For higher temperatures, liquid instead of solid product appears reducing the efficiency of the process and leading to its termination after a small number of cycles. Further to the thermodynamics, the kinetics of the reactions (1) is very important. The water splitting reaction may be actually consisted of several elementary steps. What is needed for the model development is not the exact reaction mechanism but the simplest
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
overall reaction rate expression that can describe the experimental data. An empirical water splitting kinetics scheme is proposed by Tsuji, Togawa, Wada, Sano, & Tamaura (1995) according to which the reaction rate depends only on the partial pressures of vapor and hydrogen and not on the state of the redox material. A much simpler rate expression containing a linear dependence on the vapor concentration is proposed by Chen and Hogan (2009). A fairly general reaction rate scheme has the form (Rs is the water splitting reaction rate in moles of produced hydrogen per unit area of redox material per second): Rs =
ks (1 − y)C 1 + kH C
(2)
where ks is a temperature dependent rate constant, is the maximum oxygen storage capacity of the redox material (total moles H2 produced per redox material surface area unit) and it also depends on the temperature and kH is a temperature dependent constant related to surface adsorption–desorption process. Finally, C is the vapor concentration and y is the instantaneous actual storage normalized by the maximum one and it takes values between 0 and 1. The quantity (1 − y) is proportional to the number of empty sites available for oxygen storage on the material. This rate expression is the outcome of combined theoretical considerations and experimental observations. Regarding regeneration, a constant rate has been proposed (Chen & Hogan, 2009) but the most reasonable choice is the first order dependence on the adsorbed oxygen i.e. Rr = kr y
(3)
where the reaction rate constant kr exhibits an Arrhenius temperature dependence. The water splitting and regeneration reaction constants are computed from experimental data taken in packed bed reactors with the redox material in the form of powders (Zygogianni et al., 2009). 2.2. Level II, monolith channel scale The transient operation of the monolith has in general a time scale large compared to the residence time of the gas in the monolith so a quasi-steady approximation for the flow of the gas can be considered. The one dimensional (1D) model of this scale includes the quasi-steady momentum, heat and mass balances for the gas stream (see Appendix A.1 for details). The diffusion of water vapor towards the wall has been modeled in Agrafiotis, Lorentzou, Pagkoura, Kostoglou, and Konstandopoulos (2006) using the Fick law that allows an analytical solution to the vapor mass balance in the gas phase. The use of Fick law is justified in those cases since the focus of the work was on modeling experimental campaigns with vapor molar fraction less than 10% in the feed stream. In general, considering larger fraction of vapor, the diffusion process must be described by Stefan Maxwell equations. A need arise here to develop perimeter averaged 1-dimensional models based on Stefan Maxwell equations instead on the Fick law. The key scope is to include the generality of Stefan Maxwell equations keeping the computational cost at minimum. Another pitfall of the existing 1D model is the ignorance of the radiation heat transfer mechanism. The physical properties of the gas are computed at the local cross sectional average value of the temperature but although it is a standard practice for the existing monolith applications it remains to be confirmed for the high temperature range met in the present application. Even with the above pitfalls the 1D model was able to predict the trends of the experimental results for hydrogen production cycles. An important computational aspect revealed by the 1D model is the stiffness of the set of equations regarding their integration in the time domain. Whereas in case of automotive converters the stiffness is due to the extremely high temperature dependence
1917
of the Langmuir–Hinshelwood expressions for the kinetic of the oxidation reactions (Oh & Cavendish, 1982), in the present case the stiffness originates from the heat transfer problem itself. The solar radiation enters the monolith as a flux boundary condition on the solid face of the monolith and the heat is moving slowly into the solid wall as a steep front. Unfortunately, this stiffness is inherent to the process and it is an obstacle in the use of higher level models to simulate the operation of the hydrogen operation unit. In order to examine the validity of the 1D model and to adjust it to accommodate the additional complexities of the actual problem (Stefan Maxwell mass transfer, radiation intra-channel heat transfer and high cross-sectional temperature non-uniformities), a computational fluid dynamics (CFD) oriented single channel model is developed (see Appendix A.2 for its description). The channel interior and the channel walls are discretized and the continuity equation, the momentum balances and the heat and mass conservation equations employing Stefan Maxwell diffusive flux relations, are solved in the gas phase combined with the heat conservation equation, the reaction kinetic rates and the equation for the accumulation of the oxygen in the solid phase. The computational cost of the three dimensional (3D) single channel model is very high so the effects of the 3rd dimension and of the radiation are studied carefully during a heat up of the monolith with solar radiation, up to the establishment of the steady state. Using the two dimensional (2D) equivalent problem with properly transformed variables it can be shown that the information expected from the 3D model can be also extracted using the 2D model. The CFD code admits detailed inclusion of the surface to surface radiation model. The analysis of the heat up showed that radiation acts as an additional thermal conductivity of the solid wall. This conductivity is highly depended on the temperature (monotonically increasing function of temperature). For materials with low thermal conductivity as the cordierite, the effect of the radiation on heat transfer in the monolith is important and it cannot be ignored. On the contrary the high thermal conductivity of the silicon carbide makes the contribution of the radiation negligible up to 1500 ◦ C so it can be ignored for the simulation of the cycling operation of the monolithic reactor reducing decisively the computational requirements. 2.3. Level III, entire reactor scale This is the macroscopic scale of the monolith. The need for models of this scale is due to the spatial non-uniformities of the inlet conditions or to boundary heat loses. For example, in the case of the automotive converters and diesel particulate filters (Kostoglou, Housiada, & Konstandopoulos, 2003) the non-uniformities are met mainly at the inlet flow distribution, at the inlet gas temperature distribution and they are created by the heat loses from the monolith to the environment. The above non-uniformity sources are of minor significance for the hydrogen production monolith reactor for which the primary cause of non-uniformity is the radiation strength distribution. Although, a perfectly flat heliostat would produce an image on the receiver of the size of the heliostat, for most applications each mirror segment on a heliostat is concaved slightly and canted towards the receiver. Mirror surface waviness, the gross curvature error of each heliostat, the errors associated with accurate positioning of each mirror segment and the effects of differential thermal growth tend to increase the image error. All of these add up optically to produce a higher flux density at the aim point which has a distribution pattern of a Gaussian shape typically used in simulation studies (Stine & Geyer, 2001). The detailed way to simulate the entire monolith is to discretize all the channels one by one but this approach is computationally out of question. A more reasonable approach is to simulate all channels using the 1D perimeter averaged model and allow communication between the channels through the heat transferred among their
1918
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
Fig. 1. Effect of radiation on the wall temperature profile evolution during monolith heating.
walls. Considering that a monolith typically contains a few thousands of channels this model is also computationally intractable. The usual approach for monolith scale simulations is to exploit the small value of the ratio between the channel side size and the monolith face area size and to homogenize the discrete channel set of equations to get a continuum model. The resulting homogeneous problem in which the channel level structure has been lost but it is indirectly conveyed through the appropriate transport coefficients can be solved by CFD codes as porous media. Many commercial CFD codes assume porous media to be always in thermal equilibrium between gas and solid phase. This problem is overcome by the so called dual grid approach according to which a solid phase and a porous medium are assumed to occupy the same space and connected to each other through the appropriate transport relations. A further discussion on the models of this level can be found in Appendix A.3. 3. Indicative results and discussion There are two modes of operation of the monolith: the first is the relative long heating period required for the cold monolith to reach the reaction temperature. During this period only inert gas is fed to the monolith and the only occurring phenomenon is heat transfer. The second mode is the periodic operation during which the reactions occur and the hydrogen is produced. The 3D single channel problem with radiation is too expensive computationally to be used for the simulation of the hydrogen production cycles. Comparison between the simple heating (first mode) results taken for 2D (parallel plates) and 3D (square cross section) channels revealed only quantitative differences that can be accommodated by modifying appropriately the heat transfer coefficient. So considering that the purpose of the CFD model is the adjustment of the 1D perimeter averaged single channel model, the 2D channel geometry can be employed. The inclusion of the radiation also requires significant computational effort. The wall surface temperature profile along the channel for typical heating conditions is shown in Fig. 1 in the absence and presence of radiation. The high thermal conductivity of the silicon carbide leads to a relative uniform temperature profile which does not give space to radiation for further smoothing of this profile. The effect of the radiation is negligible so radiation can be completely ignored for monoliths made by silicon carbide. The production cycles are simulated using the 2-D geometry without radiation. The simulated scenario is the following: In order to converge fast to a limit cycle, a steady state operation at 1660 K with 100% nitrogen as gas phase and all the redox material in the reduced state is assumed as initial condition. Then the solar heat flux decreases and when the temperature goes down from the 1473 K the feed changes to 80% nitrogen and 20% water vapor.
Fig. 2. Channel wall average temperature and gas outlet temperature during the periodic operation of the monolith.
When the temperature reaches the 1073 K the solar heat flux goes again to its initial value and at 1473 K the feed change again to 100% nitrogen. The convergence to a periodic steady state occurs at the second cycle. The evolution of the average wall temperature and of the gas temperature at the outlet of the monolith is presented in Fig. 2. The second cycle is the limiting one and the system has come to a periodic steady state. There is a delay between the gas and the average solid temperature. The outlet gas is in thermal equilibrium with the outflow edge of the monolith wall which always exhibits a delay with respect to the average wall temperature since the heat source is located at the inlet edge. The evolution of the concentration of the important species in the outflow is shown in Fig. 3. During the water splitting stage the vapor is completely transformed to hydrogen. When most of the redox material is in the oxidized state the rate of the reaction is reduced and a part of the vapor remains unreacted. The reaction efficiency of the water splitting process is 83% and the cycle efficiency (duration of water splitting step to total cycle duration) is 54% in the present case. The evolution of the fraction of the redox material in the reduced and oxidized states is shown in Fig. 4. The present case corresponds to an efficient operation scenario since the water splitting stage stops at oxidized state fraction y close to 0.8 and starts at y close to 0.05. The surface wall temperature profiles along the channel at several moments during the cycle are shown in Fig. 5. As it is obvious, the temperature profiles along the channel are highly non-uniform due to the combination of the strongly localized heating (entrance of the channel) and the need for periodic operation. The nonuniform conditions along the filter combined to the very strict requirements on the local temperature range (especially during the regeneration step) renders the selection and optimization of the proper operation strategy, a very demanding task. The evolution of
Fig. 3. Chemical species concentration at the outlet gas stream versus time.
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
Fig. 4. Evolution of the (complementary) fractions of oxidized (MO) and reduced (M) redox material.
1919
Fig. 7. Contour plot of the hydrogen mass fraction in the channel of the monolith at t = 720 s.
Fig. 8. Temperature (K) contour plot in the channel of the monolith at t = 720 s.
Fig. 5. Typical temperature profiles along the channel wall indicating the extent of temperature non-uniformity.
the fraction y of the oxidized redox material along the channel is shown in Fig. 6. Both oxidation and regeneration proceed from the inlet to the outlet creating highly non-uniform profiles of y. The distribution of the hydrogen in the channel during the water splitting stage is shown in Fig. 7 (at t = 720 s). The occurrence of the reaction on the wall leads to large hydrogen concentrations close to the wall. The vapor is completely consumed along the flow and the already produced hydrogen is transferred from the wall to the whole channel cross section by diffusion. The species concentration is reduced along the flow due to the gas thermal expansion as the temperature increases. The temperature contours in the channel at the same time (t = 720 s), are shown in Fig. 8. The gas is losing fast its thermal content due to its low thermal capacity and it is adjusted
Fig. 6. Evolution of profiles of the fraction y of oxidized redox material along the channel.
to follow the wall temperature. The wall temperature is uniform across the wall and varies only along the channel. The gas temperature follows the longitudal variation of the wall temperature profile retaining a self-similar transverse profile that corresponds to the asymptotic Nusselt number. The whole structure of the temperature field in Fig. 8 indicates that a 1D approach can accurately predict the heat transfer dynamics in the channel. The full scale monolith model is needed in order to study the effect of heat source non-uniformities and macroscopic geometry. A typical industrial scale setup of the monolith was studied during its heating by a heat source of Gaussian type with its maximum at the center line of the monolith and a close to zero value of flux at the circumference. Several steady state results indicating the effect of source non-uniformity on the temperature non-uniformities are presented. The temperature contours and the gas phase streamlines are shown in Fig. 9. At this temperature scale transverse nonuniformities of the temperature are not visible. The small thermal capacity of the gas in the mixing chamber result in its preheating before entering the monolith so the gas stream enters into the monolith in thermal equilibrium. The effect of the macroscopic equipment structure must be considered by channel level models that typically assume different temperature between monolith and inlet gas. The temperature profiles on three types of commercial monoliths (with 90, 180 and 270 channels per square inches – cpsi) are shown in Fig. 10. The large heat source non-uniformity has been drastically smeared-out by the high conductivity of the material leading to small temperature non-uniformities. As the size of the
Fig. 9. Temperature (K) contours (a) and gas flow streamlines (b) for a heated by non-uniform solar heat source monolith at steady state.
1920
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
Fig. 10. Steady state temperature (K) contours for monoliths (with different channel densities) heated by non-uniform solar heat source.
Fig. 12. Schematic and variables of single channel.
Appendix A. A.1. 1D single channel model
Fig. 11. Wall temperature (K) for channels at the center and at the circumference of monoliths (with different channel densities) heated by non-uniform solar heat source.
channel decreases the effective thermal conductivity of the monolith increases leading to smaller temperature non-uniformities as shown in Fig. 10. The temperature profiles along the channel for channels being at the center and at the circumference of the monolith are presented in Fig. 11. The non-uniformities at the face cross section are eliminated by the conductance along the channel as heat “diffuses” from the hotter central area to the colder circumference of the monolith. Conclusively, one could say that at least for silicon carbide the temperature non-uniformities created by the solar flux non-uniformities are small enough to be ignored in operation strategy optimization studies.
The equations that describe the single channel dynamics are presented in the following (where ˛ is the side size of the square channel, L the length of monolith and w the thickness of the wall of the monolith – see Fig. 12). The unknowns of the model are the evolved spatial (along the channel) profiles of gas phase composition and temperature, of wall temperature and of state y of the redox material. Mass conservation equations are preferably written in terms of molar flux of species instead of concentrations in order to avoid complications arising by the gas density variation (A is the redox material surface area per filter wall surface area):
A.1.1. Gas phase Species conservation: (i) H2 O (vapor) df1 = −4˛kD (c1 − c1w ) dz
(A1)
The flux to the wall must be equal to the reaction rate: kD (cH − cHw ) = ARs
(A2)
(ii) inert gas 4. Conclusions The need to optimize the operation of hydrogen production from heated monoliths process asks for advanced modeling of the process. Unfortunately, the disparate size and time scales appeared in the problem do not admit the use of a general purpose simulation tool for the demanding optimization study of the process. Different models must be used for different purposes. In this view, the three-dimensional complete monolith and housing model is necessary to study (under simple operations conditions as pure heating or steady state) the effect of heat flux or flow non-uniformity to the uniformity of the conditions at the outlet of the monolith. The detailed CFD single channel model is needed to properly adjust and validate 1D perimeter averaged single channel model under typical operation conditions. Finally, the properly adjusted 1D channel model is the appropriate scale model to be used for the extensive computations needed in view of optimization of hydrogen production unit operation. Regarding monolith material, it was shown that the high thermal conductivity of the silicon carbide simplifies considerably the modeling requirements by eliminating transverse temperature gradients and allowing ignorance of the intrachannel radiation heat transfer. The computational requirements are much higher for low conductivity materials like cordierite in which steep temperature fronts appear making necessary the consideration of transverse non-uniformity and intrachannel radiation by the model.
df2 =0 dz
(A3)
(iii) H2 df3 = 4˛ARs dz
(A4)
(iv) O2 df4 = 4˛ARr dz
(A5)
Energy conservation:
4
fi cpi
dT = 4˛kh (Tw − T ) dz
(A6)
i=1
A.1.2. Solid phase Surface coverage evolution:
dy = ARs − ARr dt
(A7)
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
Energy conservation: 2s cpw w˛
∂Tw = 2w w˛ ∂t
∂2 T ∂z 2
+ 4˛kh (T − Tw )
+ 4˛A(Hs Rs + Hr Rr )
(A8)
The concentrations in the gas phase is computed as P ci = 4 R f gT i=1 i
(A9)
In the above equations cpi is the molar specific heat capacity of the gaseous species i, cpw the specific heat capacity of the channel wall, , s the thermal conductivity and density of the wall, H are the reaction enthalpies, kh the channel to wall heat transfer coefficient and kD the channel to wall mass transfer coefficient. The molar fluxes of vapor, inert gas, hydrogen and oxygen are denoted as f1 , f2 , f3 and f4 respectively. Finally P, T denote the local pressure and temperature and Rg the gas constant. The conductivity of the porous wall is computed using typical expressions for the conductivity of composite media and the temperature dependent intrinsic conductivity of the silicon carbide. In case of siliconized silicon carbide being the monolith material, the solid phase intrinsic conductivity is computed through a weighted average of silica and silicon carbide. The mass transfer coefficient kD can be modified to account for the effect of Stefan Maxwell equations without increasing the mathematical complexity of the problem. Instead of the ternary mixture of water vapour, hydrogen and inert gas a binary system of vapour diffused in a mixture of hydrogen and inert gas is considered. Extensions of Fick’s law for binary mixtures and for diffusion in mixture of gases can be combined to give a concentration dependent diffusion coefficient which through the Sherwood number relation for square channel (asymptotic Sherwood number equal to 3) leads to the coefficient kD . The resulting relation for the effective diffusivity of vapor in the mixture of inert gas and hydrogen is Def = (xH /DwH + xin /Dwin )−1 where xH , xin are the local molar fractions of hydrogen and inert gas respectively, and DwH , Dwin are the binary diffusivities of vapor in hydrogen and inert gas respectively. A.1.3. Boundary conditions The boundary conditions for the monolith are as follows: fi (0, t) = fio (t)
(A10a)
T (0, t) = To (t)
(A10b) at z = 0
∂Tw = 0 at z = L ∂z
(A10c) (A10d)
where Q(t) is the solar power adsorbed by the monolith per channel. A.1.4. Initial condition Tw (x, 0) = Two
coupling between gas and solid phase is achieved through continuity of temperature, heat flux and reacting species fluxes. The surface to surface radiation model is considered since the gas filling the channels does not absorb radiation. Finally the temperature dependence of the gas properties is taken into account. A.3. 3D model
fi
∂Tw Q (t) − = ∂z (2w˛ + w2 )
1921
(A10e)
A.2. CFD-based detailed single channel model Only a description of the mathematical problem for this level model will be given here to avoid copying of equations described in the manual of a commercial CFD code (specifically the code FLUENT has been used here). The governing equations are the continuity equation, the momentum conservation equations, the Stefan Maxwell relations for the diffusive fluxes, the gas phase energy conservation equation and the solid phase thermal conservation equation. The surface reaction feature of the code is used for the consideration of water splitting and regeneration reactions. The
The three dimensional model results from a homogenization procedure during which the geometrical details of the monolith are lost and the solid and gas phases are considered as interpenetrating continua. The governing equations are the same with those of the 1dimensional single channel model but with two differences: (i) Two additional spatial independent variables are needed (x, y) rendering all the depended variables functions of x, y, z, t. (ii) A mesoscopic conduction term has to been added. The effective conductivity tensor has two components one for the axial direction and one for the transverse direction. Derivation of these components for monoliths with square channels can be found elsewhere (Groppi & Tronconi, 1996; Konstandopoulos, Kostoglou, Vlachos, & Kladopoulou, 2007). The resulting system of equations can be solved by assuming representative channels on the monolith cross-section, solving them as 1D problems, find the heat sources for each channel and update the solid phase heat transfer equation with these sources. Alternatively, the CFD code FLUENT can be used through the so called dual grid approach considering a porous medium and a compact solid properly connected to each other. The reactor scale results presented here were taken employed this approach. References Agrafiotis, C., Roeb, M., Konstandopoulos, A. G., Nalbandian, L., Zaspalis, V. T., & Sattler, C. (2005). Solar water splitting for hydrogen production with monolithic reactors. Journal of Solar Energy Engineering, 79, 409–421. Agrafiotis, C., Lorentzou, S., Pagkoura, C., Kostoglou, M., & Konstandopoulos, A. G. (2006). Advanced monolithic reactors for hydrogen generation from solar water splitting. In Seville: Proceedings of solarPACES 13th international symposium on concentrated solar power and chemical energy technologies 2006 Seville, Spain, June 20–23. Agrafiotis, C., Pagkoura, C., Lorenzou, C., Kostoglou, M., & Konstandopoulos, A. G. (2007). Hydrogen production in solar reactors. Catalysis Today, 127, 265–277. Allendorf, M. D., Diver, R. B., Siegel, N. P., & Miller, J. E. (2008). Two step water splitting using mixed-metal ferrites: Thermodynamic analysis and characterization of synthesized materials. Energy & Fuels, 22, 4115–4124. Chen, K. S., & Hogan, R. E. (2009). A two phase model for solar thermochemical water splitting with FeO/Fe3 O4 . In Proceedings of ES2009. 2009 ASME 3rd international conference of energy sustainability San Francisco, California, July 19–23. Deen, W. M. (1998). Analysis of transport phenomena. Oxford University Press: New York. Dersch, J., Mathijssen, A., Roeb, M., & Sattler, C. (2006). Modelling of a solar thermal reactor for hydrogen generation. Modelica. Dubien, C., Schweich, D., Mabilon, G., Martin, B., & Prigent, M. (1998). Three-way catalytic converter modelling: Fast- and slow-oxidizing hydrocarbons, inhibiting species, and steam-reforming reaction. Chemical Engineering Science, 53, 471–481. Groppi, G., & Tronconi, E. (1996). Continuous vs. discrete models of non-adiabatic monolith catalysts. AIChE Journal, 42, 2382–2387. Koci, P., Kubicek, M., & Marek, M. (2004). Modeling of three-way-catalyst monolith converters with microkinetics and diffusion in the washcoat. Industrial & Engineering Chemistry Research, 43, 4503–4510. Kodama, T., & Gokon, N. (2007). Thermochemical cycles for high-temperature solar hydrogen production. Chemical Reviews, 107, 4048–4077. Konstandopoulos, A. G., Kostoglou, M., Vlachos, N., & Kladopoulou, E. (2007). Advances in the science and technology of diesel particulate filter simulation. In G. B. Marin (Ed.), Advances in chemical engineering, Special issue on automotive emission control (pp. 213–275). Academic Press. Kostoglou, M., & Konstandopoulos, A. G. (2005). Effect of soot layer microstructure on diesel particulate filter regeneration. AIChE Journal, 51, 2534–2546. Kostoglou, M., Housiada, P., & Konstandopoulos, A. G. (2003). Multi-channel simulation of regeneration in honeycomb monolithic diesel particulate filters. Chemical Engineering Science, 58, 3273–3283. Mladenov, N., Koop, J., Tischer, S., & Deutschmann, O. (2010). Modeling of transport and chemistry in channel flows of automotive catalytic converters. Chemical Engineering Science, 65, 812–826. Oh, S. H., & Cavendish, J. C. (1982). Transients of monolithic catalytic converters: Response to step changes in feedstream temperature as related to controlling
1922
M. Kostoglou et al. / Computers and Chemical Engineering 35 (2011) 1915–1922
automobile emissions. Industrial & Engineering Chemistry Product Research and Development, 21, 29–37. Roeb, M., Sattler, C., Klüser, R., Monnerie, N., de Oliveira, L., & Konstandopoulos, A. G. (2006). Solar hydrogen production by a two-step cycle based on mixed iron oxides. Journal of Solar Energy Engineering, 128, 125–133. Roeb, M., Neises, M., Säck, J. P., Rietbrock, P., Monnerie, N., & Dersch, J. (2008). Operational strategy of a two-step thermochemical process for solar hydrogen production. International Journal of Hydrogen Energy, doi:10.1016/j.ijhydene.2008.08.049 Roeb, M., Gathmann, N., Neises, M., Sattler, C., & Pitz-Paal, R. (2009). Thermodynamic analysis of the two-step solar water splitting with mixes iron oxides. International Journal of Energy Research, doi:10.1002/er.1513 Siemund, S., Leclerc, J. P., Schweich, D., Prigent, M., & Castagna, F. (1996). Three-way monolithic converter: Simulations versus experiments. Chemical Engineering Science, 51, 3709–3720. Steinfeld, A., Sanders, S., & Palumbo, R. (1999). Design aspects of solar thermochemical engineering – A case study: Two-step water-splitting cycle using the Fe3 O4 /FeO redox system. Solar Energy, 65, 43–53.
Stine, W. B., & Geyer, M. (2001). Power from the sun (electronic edition). Tischer, S., & Deutschmann, O. (2005). Recent advances in numerical modeling of catalytic monolith reactors. Catalysis Today, 105, 407– 413. Tischer, S., Correa, C., & Deutschmann, O. (2001). Transient three-dimensional simulations of a catalytic combustion monolith using detailed models for heterogeneous and homogeneous reactions and transport phenomena. Catalysis Today, 69, 57–62. Tsuji, M., Togawa, T., Wada, Y., Sano, T., & Tamaura, Y. (1995). Kinetic study of the formation of cation-excess magnetite. Journal of the Chemical Society – Faraday Transactions, 91, 1533–1538. Young, L. C., & Finlayson, B. A. (1976). Mathematical-models of monolith catalytic-converter. Application to automobile exhaust. AIChE Journal, 22, 343– 353. Zygogianni, A., Lorentzou, S., Pagkoura, C., Lekkos, C., Kostoglou, M., & Agrafiotis, C. (2009). On the modeling of solar monolith reactors for hydrogen production. In Proceedings of SolarPACES 2009 conference Berlin.