A one-dimensional multicomponent two-fluid model of a reacting packed bed including mass, momentum and energy interphase transfer

International Journal of Multiphase Flow 57 (2013) 10–28

Contents lists available at SciVerse ScienceDirect

International Journal of Multiphase Flow j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j m u l fl o w

A one-dimensional multicomponent two-fluid model of a reacting packed bed including mass, momentum and energy interphase transfer Robert-Jan Koopmans a, John S. Shrimpton a,⇑, Graham T. Roberts a, Antony J. Musker a,b a b

Aerodynamics and Flight Mechanics Group, University of Southampton, 1 University Road, Southampton SO17 1BJ, United Kingdom DELTACAT Ltd., Kintyre House, 70 High Street, Fareham PO16 7BB, United Kingdom

a r t i c l e

i n f o

Article history: Received 25 February 2013 Received in revised form 13 June 2013 Accepted 14 June 2013 Available online 27 June 2013 Keywords: Two-fluid flow Multicomponent Chemical reaction Evaporation Packed bed Hydrogen peroxide

a b s t r a c t In this paper a multicomponent two-phase flow numerical model for reacting packed beds is presented, where one phase is compressible. The model is applied to a catalyst bed which decomposes highly concentrated hydrogen peroxide into oxygen and water vapour. Using dimensional analysis a number of simplifications to the governing equations are made. A pressure based finite volume approach is used to solve the resulting equations where the pressure–velocity–density coupling is established with a SIMPLE algorithm for two-fluid flows. The model is validated against experimental data and shows in general a reasonable agreement given the complexity of the physics. It is shown that with the two-fluid formulation characteristics are predicted that cannot be reproduced with mixture models. Simulations reveal that the largest part of the total mass transfer takes place in a very small section of the catalyst bed at the point where the liquid reaches the boiling point and that most of this transfer is due to evaporation. It is further shown that a thermal description of the packing of the catalyst bed would improve the prediction of temperature distribution in the two-phase regime. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Hydrogen peroxide was extensively used as rocket propellant in the 1950s and 1960s, but was gradually replaced after by propellants with better performance characteristics, such as hydrazine and combinations of MMH (CH6N2) and UDMH (C2H8N2) with NTO (N2O4) (Wernimont et al., 1999). Recently the interest in the use of hydrogen peroxide as fuel for rocket motors has revived (Kappenstein, 2008) driven by the toxic and carcinogenic properties of rocket fuels currently used and the associated potential cost savings by using less harmful propellants. A hydrogen peroxide based rocket motor relies on the decomposition of the peroxide in a decomposition chamber where it is brought into contact with a catalyst. The catalyst can be in liquid form, in which case it is injected separately from peroxide into the decomposition chamber (Musker et al., 2004), or in solid form as silver screens (Helms et al., 2002), pellets (Pasini et al., 2008) or monoliths (Fanelli et al., 2011). The exothermic decomposition into oxygen and water vapour causes the temperature to rise, which leads to increased reaction and evaporation rates. The gas that leaves the decomposition chamber can reach a temperature of nearly 1270 K for anhydrous peroxide at standard atmospheric ⇑ Corresponding author. Tel.: +44 (0)23 8059 4894; fax: +44 (0)23 8059 3058. E-mail addresses: [email protected] (J.S. Shrimpton), [email protected] (G.T. Roberts). 0301-9322/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijmultiphaseflow.2013.06.005

conditions. The hot gases can be either exhausted through a nozzle, for monopropellant rockets, or mixed with a fuel for combustion, for bipropellant rockets. Despite the renewed interest in hydrogen peroxide the number of flow models available in the literature describing the processes in the decomposition chamber for porous solid catalysts is limited, probably due to the severe complexity of the system. A simple model was presented by Johnson et al. (2000), where the decomposition chamber was modelled by two control volumes: the control volume adjacent to the injector containing all the liquid and a control volume downstream containing the gas portion. The purpose was to investigate pressure oscillation experienced during testing. Liquid to gas mass transfer was assumed to have taken place after the liquid was in the decomposition chamber for a fixed amount of time. No distinction was made between decomposition and evaporation. A one-dimensional model was developed by Zhou and Hitt (2003). In this model the gas and liquid phase were modelled as a mixture. Contrary to Johnson et al. (2000) they considered decomposition and evaporation separately which were driven by the mixture temperature instead of a fixed residence time of the liquid. They only considered catalytic decomposition, which was modelled with Arrhenius kinetics. For the evaporation of the liquid they used a staged approach in which water and peroxide evaporated at their respective normal boiling points. Their results showed that this unphysical assumption only affected a very small

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

fraction of the catalyst bed at the inlet. Validation against experimental data was not given. Their model focussed on a micropropulsion system with a catalyst bed length in the order of hundreds of micrometers. The catalyst bed itself had the form of a monolith. They performed parametric studies in which they varied the Arrhenius parameters and heat loss to the environment to investigate the influence on the required bed length to ensure complete decomposition. Bonifacio and Sorge (2006) investigated the transient behaviour of a monolithic catalyst bed and the thermal response of the walls in particular. For this purpose they developed a lumped-parameter model in which catalytic decomposition took place at the wall. Heat loss to the environment was not taken into account. To simplify the analysis they only modelled the gas phase and ignored the liquid phase. They found out that during the transient startup the wall could reach higher temperatures than the adiabatic decomposition temperature before a steady state temperature was reached equal to the adiabatic decomposition temperature. This effect was dependent on the relative importance of mass and heat transfer by diffusion. In a model developed by Corpening et al. (2006) special attention was paid to thermal decomposition in the gas phase and evaporation of highly concentrated peroxide droplets. Droplets were assumed to be surrounded by a hot gas with a temperature of 1200 °C and higher. The evaporation of these droplets was modelled with the D2 law under the assumption that the vapour coming from the droplet surface had the same composition as the liquid. They argued that this assumption was allowed as the droplets were small, the surrounding gas temperature high and therefore the evaporation process was very fast. Pasini et al. (2010) developed a one-dimensional mixture model for a pellet based catalyst bed. Contrary to the work by previous researchers they split the decomposition process into adsorption onto and desorption from the catalyst pellets. They also took into account the difference in peroxide concentration between the bulk flow and the fluid near the pellet surface by means of the Reynolds analogy. Similar to Corpening et al. (2006) they assumed for evaporation that the vapour coming from the droplet surface had the same composition as the liquid. They also explicitly modelled the viscous interaction between the fluid and the catalyst pellets by means of the Ergun equation. Two-dimensional models were mentioned by Zhou and Hitt (2005) and Bonifacio and Sorge (2006). Zhou and Hitt compared their earlier developed one-dimensional model (Zhou and Hitt, 2003) with a two-dimensional CFD model. They reported a large difference in required bed length to ensure complete decomposition between their one-dimensional model and the CFD model. Bonifacio and Russo Sorge compared the results from their lumped-parameter model with a two-dimensional CFD simulation performed in Fluent. They found a good agreement between the two models. They also mentioned that they had problems in reaching convergence with the CFD model. Unfortunately, in both cases no further details were provided about the CFD model. Employing simple mixture models to describe the processes inside the catalyst bed, such as has been done in all previous work discussed above, implicitly assumes equilibrium between the phases (Ishii and Hibiki, 2006). However, there is reason to believe that such equilibrium does not exist. For example, the adiabatic decomposition temperature of 80 wt.% hydrogen peroxide initially at room temperature is about 780 K, while the critical temperature for peroxide and water is 739.5 K and 647.3 K respectively (Manatt and Manatt, 2004). In this paper a one-dimensional two-fluid model is presented describing the decomposition and evaporation of hydrogen peroxide in a pellet-based catalyst bed in which the pellets are coated with the catalytic active phase. The aim of this model is to

11

investigate global characteristics or first-order approximations of a decomposing flow of highly concentrated hydrogen peroxide. Both fluid phases are multi-component, the liquid incompressible and the gas phase compressible. For this purpose the Navier– Stokes equations are rewritten in such a way that the model can be solved for volume fractions from near 0 to near 1. Special attention is paid to the development of the source terms describing decomposition, evaporation and the interaction with the catalyst pellets. The resulting simulations are compared with experimental data. Although the model is developed for hydrogen peroxide based rocket applications, the model is generally applicable to two-fluid problems in packed beds with interphase mass, momentum and heat transfer. Examples include hydrogenation in the food industry for the production of margarine from vegetable oils and the Fischer–Tropsch process for the production of synthetic fuels. This paper starts with a presentation of the two-fluid equations in the most general form. Special attention is being paid to the formulation of the mass, momentum and energy source terms. Then a non-dimensional analysis is performed from which a set of simplified equations is derived. Special numerical considerations are discussed that permit stable solution at realistic reaction rates. The section after that discusses the experimental validation of the model. Finally, some results of a typical simulation are discussed where the benefits of a two-fluid model over mixture models will be emphasised. 2. Model formulation In this section the conservation equations for the one-dimensional flow are presented. As the system under consideration is highly complex, assumptions and simplifications have to be made to make the problem tractable. To justify these, first the complete set of conservation equations are given after which for each term the relative importance will be assessed. This will result in a simplified set of equations that are solved. 2.1. Governing equations Given that the aim of the current model is to investigate firstorder effects, variations in longitudinal direction only are considered: it is assumed there is no variation of properties in the radial direction. Consequently, the average over the cross section is equal to the local conditions at any radial position. This implies that there are no gradients in radial direction for all variables. The density, q, velocity, u, static enthalpy, h, and species mass fraction, Y, for each component a can then be determined by solving the continuity, momentum, enthalpy and species equations for each phase. The resulting set of instantaneous equations for each phase k is given by Ishii and Hibiki (2006), where z is in axial direction and i refers to an interface quantity

@ k qk @ k qk uk þ ¼ Ci;k ; @t @z @ k qk uk @ k qk u2k @p @ k sk;zz @ k þ ¼ k k þ þ ðpi  pk Þ @t @z @z @z @z @ k þ ui;k Ci;k  si;zz þ M i;k ; @z   @ k qk hk @ k qk uk hk @ @T Dp @u þ ¼ k kk k þ  k þ k sk;zz k @z @t @z @z Dt @z Dk @ k  si;zz ðui;k  uk Þ þ ðpk  pi Þ Dt @z þ hi;k Ci;k þ ui;k M i;k þ Hi;k ;

ð1Þ

ð2Þ

ð3Þ

12

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

@ k qk Y k;a @ k qk uk Y k;a @ þ ¼ @z @t @z



k Dk;a

 @Y k;a þ Sk;a : @z

ð4Þ

Here p is the pressure, szz the normal stress,  the fluid volume fraction ranging from 0 to 1, T the temperature, k thermal conductivity and Da the binary diffusion coefficient. The source terms are given by Ci,k, representing the mass source for each phase; Mi,k, a momentum source due to interaction of the phases; Hi,k, describing the energy transfer due to heat transfer between the phases, and Sk,a, the mass source for each species a. Note that summation of the mass P sources over all the phases k gives k Ci;k ¼ 0 The mass source itself can be split into a contribution from decomposition and evaporation, Ci;k ¼ Ci;k dec þ Ci;k ev ap . The energy transfer due to interface mass transfer is described by the 6th term on the RHS of Eq. (3). Eq. (4) accounts for the multiple species in each phase. Here the liquid phase consists of two species, hydrogen peroxide and water, and the gas phase of four species, hydrogen peroxide vapour, water vapour, oxygen and air. The reason to have air included will become apparent below. The presence of the pellets in the catalyst bed manifests itself via the momentum source term in Eq. (2), which will be further discussed in Section 2.2.3, and via the reduced amount of volume available for fluid flow. The void fraction f is dependent on the catalyst bed diameter and the pellet geometry and dimensions. Normally, the presence of the wall leads to an ordering effect near the wall (Bey and Eigenberger, 1997). For catalyst beds with a smaller diameter or for larger pellets, these effects extend further towards the center of the catalyst bed. Bey and Eigenberger (1997) presented an equation for the average void fraction as a function of the bed diameter and pellet diameter. For cylindrical pellets it is given by

f

¼ 1  s

 2 dps dps ¼ 0:36 þ 0:1 þ 0:7 ; dt dt

ð5Þ

where s is the solid volume fraction in the bed, dt is the diameter of the catalyst bed and dps the diameter of a sphere with an equivalent volume of a catalyst pellet. It is valid for dps/dt 6 0.6. Any heat transfer between the catalyst pellets and the fluids as well as the thermal mass of the pellets is not accounted for in the present model. The implications of this assumption are discussed in the validation. 2.2. Source terms In the current model four sources can be distinguished: mass sources originating from decomposition of hydrogen peroxide and from evaporation of the liquid, a momentum source accounting for the friction between the different phases in the catalyst bed and heat transfer between the fluids due to a temperature difference. Note that part of the liquid decomposition results in interfacial mass transfer, represented by Ci in Eq. (1), and part only affects the species concentration as no phase change occurs, represented by Sa in Eq. (4). The energy that is released due to decomposition and the energy consumed due to evaporation are driven by mass transfer and will therefore not be discussed separately. 2.2.1. Decomposition Decomposition of hydrogen peroxide is generally given by (Mol et al., 2008)

1 H2 O2 ! H2 O þ O2 ; 2

ð6Þ

which is highly exothermic. Decomposition can take place in two different ways: catalytically (Chou and Huang, 1999), in which case the reaction takes place at, and is strongly promoted by, the catalyst surface; and thermally, in which case the temperature

determines the rate of decomposition (Pearson et al., 2003). To model catalytic decomposition external and internal diffusion as well as adsorption and desorption should be accounted for. Oehmichen et al. (2010) showed that for the liquid phase at high peroxide concentrations external as well as internal diffusion can be neglected. Although not investigated, the same will be assumed for the gas phase. The catalytic decomposition process can now be written as ka

kreac

A þ site  A  site ! products; kd

ð7Þ

where ka and kd are the adsorption and desorption rate constants, respectively, and kreac the reaction rate constant for decomposition. The adsorption and desorption processes are assumed to be very fast and therefore can be considered to be in equilibrium. In that case it can be described by the Langmuir isotherm (Atkins and De Paula, 2006)

ka ½H2 O2 Nð1  hÞ  kd Nh;

ð8Þ

where [H2O2] is the concentration of hydrogen peroxide, N the number of active sites on the catalyst surface and h the relative occupancy of N. Rewriting this as an expression for h results in

h¼

K r ½H2 O2  ; 1 þ K r ½H2 O2 

ð9Þ

where Kr is the ratio of the adsorption and desorption constants, Kr = ka/kd. If it is assumed that Kreac  kde, as according Pasini et al. (2010) is allowed, then the reaction rate r_ can be expressed as

r_ S ¼ kreac N

K r ½H2 O2  : 1 þ K r ½H2 O2 

ð10Þ

The subscript S refers to the conditions at the catalyst surface and the minus sign accounts for the fact that hydrogen peroxide is being consumed. If it is further assumed that decomposition rate constant can be written as an Arrhenius equation (Pasini et al., 2010) and that the equilibrium constant is very small, then the reaction rate can be written as EA

r_ S ¼ A0 e Rc T NK r ½H2 O2 :

ð11Þ

A0 is the pre-exponential factor, EA the activation energy, Rc the universal gas constant and T the temperature of the phase for which the desorption rate constant is being determined. Eqs. (8) and (9) were originally developed for gas phase reactions and show large deviations from experimental data when applied to liquids. Several researchers have proposed different modifications to Eq. (9) for liquids by introducing more empirical constants and factors. As these are not known a priori the equilibrium constant Kr will be the unknown variable for which a value has to determined experimentally, see Section 4. Another problem is that during experimental determination of the Arrhenius parameters the number of active sites N is not accounted for. Consequently, in practice the value for A0 includes the number of active sites. The same approach will be adopted here. The volumetric reaction rate r_ can now be expressed as EA

r_ H2 O2 ;S ¼ A0 e Rc T ½H2 O2  Asl K r ;

ð12Þ

where Asl is the solid–liquid interfacial surface area per unit of volume, enabling the conversion from a surface reaction rate to a volumetric reaction rate. The model for the interfacial area is described in Appendix B. Thermal decomposition on the other hand is a volumetric reaction of the reaction rate is given by EA

r_ H2 O2 ;V ¼ A0 e Rc T ½H2 O2 ;

ð13Þ

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

Multiplying Eqs. (12) and (13) by the molar mass of peroxide results in a reaction rate in kg m3 s1. The source term Ci,kdec can now be written as

Ci;k dec ¼ ðr_ H2 O2 ;S þ r_ H2 O2 ;V ÞM;

ð14Þ

where M is the molar mass. 2.2.2. Evaporation The boiling point of a multicomponent liquid is reached whenever the sum of the vapour pressure of each liquid component is equal to or higher than the surrounding pressure. The vapour pressure of a component is a function of the activity coefficient c, mole fraction in the liquid phase x and the vapour pressure of the pure component p° and for component a can be written as (Hatsopoulos and Keenan, 1981)

pa;v ap ¼ ca xa pa :

ð15Þ

Relations for the activity coefficient and vapour pressure of the pure components water and hydrogen peroxide were presented by Manatt and Manatt (2004). For a given pressure and known vapour pressure of a component the corresponding mole fraction y in the gas phase, under assumption of equilibrium between the phases, can be found with Dalton’s law

ya ¼

pa;v ap : p

ð16Þ

The difference in mole fraction between the liquid and the gas phase drives the evaporation of the liquid. The multicomponent evaporation model of Brenn et al. (2007) is used and is listed in full in Appendix C. Below several adjustments to the original model are discussed. The model which was presented by Brenn et al. (2007) is an extension of the model for multicomponent liquids and originally presented by Abramzon and Sirignano (1989). The model developed by Brenn et al. (2007) determines the evaporation rate for each component as if the liquid consists of just a single component. The rates are then multiplied by a volume equivalent partial radius, rV,a, which corresponds to the volume fraction of the component for which the evaporation rate is calculated. Brenn et al. (2007) found a relation of r V;a ¼ 0:5dl w1=3 a to give a good agreement with experimental data, where dl is the droplet diameter and wa the volume fraction of component a in the liquid mixture. It is a scaling factor for the surface area of the droplet that can be interpreted as the reduction in surface area coverage by component a due to the presence of other components at the surface of the droplet. The total rate of evaporation is then the sum of the evaporation rates of the individual components and is expressed as

_ ev ap ¼ m

X X _ a;ev ap ¼ 2p r V;a qk Da Sha lnð1 þ BM;a Þ: m a

ð17Þ

a

D is the binary diffusion coefficient, Sh⁄ the modified Sherwoord number and BM;a is the Spalding mass transfer number for component a and is defined as

BM;a ¼

Y a;s  Y a;1 P ; 1  a Y a;s

ð18Þ

where Ya,s is the gas mass fraction of component a at the droplet surface and Ya,1 the mass fraction of component a far away from the droplet surface. When the liquid mixture approaches the boiling P point a Y a;s ! 1 and consequently BM;a ! 1 and thus the energy released by decomposition is consumed by evaporation. Note that in the original model of Brenn et al. (2007) the denominator is defined as 1  Ya,s. As a result evaporation will not be enough to consume all the energy and the liquid temperature will increase to unphysically high values. Because Brenn et al. (2007) validated their

13

model only for temperatures well below the boiling point, this problem was not apparent. Eq. (17) describes the evaporation of a single droplet. To convert _ ev ap to an evaporation rate per unit volume the droplet diameter m dl should be chosen such that the droplet surface area is equal to the gas–liquid interfacial surface area as determined by the interfacial area model, see Appendix B. This means that the unit for _ ev ap in Eq. (17) is in kg s1 per unit cell volume. Multiplying this m  , gives the correct unit for the evapby the number of cell per m3, n orative mass source term

_ k;ev ap n : Ci;k ev ap ¼ m

ð19Þ

It should be noted that in reality a mix of film boiling and pool boiling and evaporation from the hot liquid surface is taking place. The combined effect will result in different local evaporation rates. To account for this the volume equivalent radius in Eq. (17) will be replaced by a constant, for which the value will be determined experimentally, multiplied by the droplet diameter; see Section 4 below. For the heat transfer between the liquid and the gas the nonequilibrium Langmuir–Knudsen model is used. In a comparison study of different droplet evaporation models by Miller et al. (1998) it was shown that this model gives predictions that are much closer to experimental values than heat transfer models based on equilibrium between the phases. The rate of change of energy of the droplet, h_ l , is formulated as

_ ev ap ; h_ l ¼ Gkg Nu ðT g  T l Þ þ DHlg m

ð20Þ

where the first terms describes the heat transfer between the phases and the second term the latent heat effect. Miller et al. (1998) showed that the improved predictions are due to the variable G which accounts for the effect of the mass transfer on the heat transfer. They also showed that when using this factor in the original model by Abramzon and Sirignano, the predictions are significantly improved. G accounts for the heat transfer reduction due to evaporation and is the analytical solution derived by Bellan and Summerfield (1978) of the Langmuir–Knudsen evaporation law. It is defined as

G¼

b ; eb  1

ð21Þ

where b is defined as

b¼

0:5cp _ ev ap : m pdl kg

ð22Þ

Here cp is the specific heat at constant pressure of the surrounding gas. Note that the driving force of first term on the RHS of Eq. (20) is the temperature difference between the fluid phases, while the driving force of the last term of Eq. (20) is the evaporation rate. For the source terms in Eq. (3) this means that

hi;k ¼ DHlg ;

ð23Þ

; Hi;k ¼ Gkg Nu ðT g  T l Þn

ð24Þ

 in Eq. (24) is used for the same reason as in Eq. (19). where n 2.2.3. Interfacial friction The interfacial friction is normally determined as a function of the relative velocity between the phases, the interfacial area and the drag coefficient (Ishii and Hibiki, 2006). The interfacial area and the drag coefficient are dependent on the volume fraction and the type of flow. In the system under consideration the gas volume fraction changes from zero to one and consequently the flow type and drag coefficient changes. The situation is further compli-

14

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

cated by the presence of catalyst pellets in the catalyst bed. For this reason an empirical model will be employed. The model computes the pressure drop per unit length and is based on the Lockhart–Martinelli parameter for two-phase flows in pipes. Zeigarnik and Kalmykov (1985) developed from this a pressure drop model for two-phase boiling flow in pebble beds. Sorokin (2008) generalised it to two-phase flows in pebble beds and compared numerical results with experimental data for mass fluxes up to 27 kg m2 s1. This showed a good agreement between experimental data and numerical results for mass fluxes up to 17 kg m2 s1 and a slight overestimation for higher mass fluxes. Typical mass fluxes for the current application are from about 50 to 300 kg m2 s1 and higher (Palmer et al., 2011a). Zeigarnik and Kalmykov (1985) and Sorokin (2008) used the Ergun equation (Bird et al., 2002) as the basis of their model. However, as the mass fluxes are expected to be a factor 2 or more higher than the range for which Sorokin validated his model, instead of the Ergun equation, which is suited for a range of 101 < Re < 103, the Tallmadge equation will be used. It is an extension of the Ergun equation for higher Reynolds numbers, suited for 101 < Re < 105 (Tallmadge, 1970), and is given by

  Dp 150 4:2 1  f qu20 ¼ þ 1=6 ; L Re 3f dp Re

ð25Þ

where Re is defined as Re = qu0dp/(l (1  f)) where l is the dynamic viscosity,  the void fraction and u the velocity. The form of the Reynolds number will be further discussed in the next subsection. Note, that in the original Tallmadge equation, and for that matter in the Ergun equation, the superficial velocity is used. The true velocity used in Eq. (25) is related to the superficial velocity by u0 = fu (Bird et al., 2002), where u0 is the superficial velocity. The complete pressure drop model for a two-phase flow in a packed bed is given in Appendix D. The momentum source term Mi,k in Eq. (2) is now be written as

M i;k ¼

Dp : L

ð26Þ

2.3. Non-dimensional analysis As the current model will be used for first-order approximations only it will be assumed that both phases share the same pressure field (so pi  pk = 0). To simplify the problem further terms that have only a minor influence will be omitted. For this purpose Eqs. (1)–(4) are first made non-dimensional and subsequently an order of magnitude analysis is carried out. Different geometrical scales can be identified for the system under consideration. The catalyst bed used for experimental validation has a diameter of 0.016 m and a length of 0.096 m. Typical dimensions of the cylindrical catalyst pellets are about 0.003 m for both diameter and length. The geometrical scaling parameter should be chosen such that geometry related quantities scale proportionally. This is the case for the void fraction f, Eq. (5), the pressure drop over the bed, Eqs. (25) and (D), and the interfacial area, Appendix B. Based on these equations dp/(1  f) is chosen as an appropriate scaling factor. Table 1 shows how the parameters are non-dimensionalised, where the subscript 0 refers to reference conditions and the superscript ⁄ to the non-dimensional form of the parameter. Note that the volume fraction  and mass fractions Y are already non-dimensional and therefore does not require further treatment. The non-dimensional continuity equation can now be written as

@ k qk @ k qk uk þ ¼ DaI Ci;k ; @t  @z

ð27Þ

Table 1 Parameter normalisation. Parameter

Unit

Normalised

q ¼ qq0

3

q

kg m

u

m s1

u ¼ uu0

z

m

z ¼

dp

m

dp ¼

t

s

p (incompressible)

N m2



zð1f Þ dp;0 dp ð1f Þ dp;0

tu ð1 Þ t ¼ 0 dp;0 f  p ¼ q pu2 0 0 p ¼ pp 0  h ¼ hh0  k ¼ kk0 c cp ¼ cp;0p 

p (compressible)

N m2

h

J kg1

k

W m1 K1

cp

J kg1 K1

l s

N s m2 N m2

l ¼ ll0 s ¼ ss0

D

m2 s1

D ¼ DD0

DaI is the first Damköhler number and is the ratio of the mass transfer rate due to decomposition and evaporation and the mass transfer rate by advection. Strictly speaking decomposition is occurring catalytically as well as thermally. The first type is a surface reaction while the second one is a volumetric reaction. In the case of surface reactions it is more appropriate to use the second Damköhler number defined as the ratio of the mass transfer rate due to decomposition and evaporation and the mass transfer rate by diffusion. However, as the dimension of the source term is in kg m3 s1, any surface reaction has to be converted into a volumetric reaction. The non-dimensional form of the momentum equation is given by 2 @ k qk uk @ k qk uk @p 1 @ k sk;zz 1 þ ¼ Ak  þ    Re @z Re @t @z @z þ DaI ui;k Ci;k þ M i;k :

si;zz

@ k @z ð28Þ

The factor A is either 1, when the incompressible phase is considered, or Eu, when the compressible phase is considered, where Eu is the Euler number, which is the ratio of the static pressure to the dynamic pressure, p/(qu2). The reason for this difference is the way in which the pressure is non-dimensionalised as is shown in Table 1. Re is the Reynolds number, defined as Re = qudp/(l(1  f)). The non-dimensional enthalpy equation can be written as 



@ k qk hk @ k qk uk hk 1 @ Pr þ ¼ Re @z @t  @z @u Dp Ec þB  þ k sk;zz k Re Dt @z  þ DaI hi;k Ci;k þ Ec ui;k M i;k 

Ec Re

si;zz ui;k



k kk

@T k @z



ð29Þ

@ k þ DaIII Hi;k : @z

Here Pr is the Prandtl number, given by cpl/k, and Ec the Eckert number, which is a ratio of the kinetic energy to the thermal energy u2/(cpT). DaIII is the third Damköhler number and is equivalent to DaI except that energy is considered instead of mass: the ratio of the rate of heat generation due to decomposition and evaporation and the heat transfer rate by advection. The factor B is either Ec in case of incompressible flow or (c  1)Ma2 in case of compressible flow. The reason is the same as explained above for the momentum equation. The non-dimensional version of the species equation is given by

@ k qk Y k;a @ k qk uk Y k;a 1 @ þ ¼ Re Sc @z @t @z



k Dk;a

 @Y k;a þ DaI Sk;a ; @z

where Sc is the Schmidt number, defined as

l=ðqDÞ.

ð30Þ

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

Table 2 shows the result of the order of magnitude analysis of the non-dimensional numbers. The first column gives an overview of the orders of magnitude for each variable. For some variables a range is given as the order changes with position in the catalyst bed. The second column shows the maximum order of magnitude for each non-dimensional number. Applying the order of magnitude numbers to the non-dimensional conservation equations shows that for the momentum equation the normal stresses can be neglected as these terms are at least 4 orders of magnitude smaller than the other terms on the RHS of Eq. (28). It also shows that for the gas phase the first term on the RHS is the dominating term. As far as the enthalpy equation is concerned, all the terms containing viscous stresses can be safely ignored as the order of magnitude of Ec/Re is lower than 8. For the same reason the heat conduction term will be ignored as 1/ RePr is of order 4. Also the work due to interfacial drag term, ui;k M i;k , will be ignored as its overall impact is small. According to the same reasoning also the third term on the RHS could be neglected as Ma2 = O(4) and Ecl = O(6). However, the term Dp⁄/ Dt⁄ can become very large during the first part of the transient simulation and for small time steps. This term will therefore be included in the simulations. Finally, the diffusion term in the species equation can be neglected as well, since the second term is much larger. Due to the reduced number of terms the momentum, enthalpy and species conservation equations become much simpler. They can be written as respectively

@ k qk uk @ k q þ @t @z

2 k uk

ð31Þ

@ k qk hk @ k qk uk hk Dp þ ¼ k þ hi;k Ci;k þ Hi;k ; Dt @t @z

ð32Þ

@ k qk Y k;a @ k qk uk Y k;a þ ¼ Sk;a : @t @z

ð33Þ

2.4. Special numerical considerations To ensure that the problem is well-posed two situations require special attention: the behaviour of the equations when one of the volume fractions tends to zero and the bounding of the mass and

Table 2 Order of magnitude of non-dimensional numbers. Variable

Max. magnitude

cp = O(3) dp = O(3) Da;g ¼ Oð5Þ Da;l ¼ Oð9Þ kg = O(2) kl = O(0) 000 kg ¼ Oð5Þ

DaI = O(2) DaIII = O(2) Ecg = O(4) Ecl = O(6) Eug = O(4) Ma = O(2) Reg = O(4)

000

kl ¼ Oð1Þ p = O(6) R = O(2) Tg = O(2)  O(3) Tl = O(2) ug = O(1)  O(1) ul = O(1)  O(0) z = O(1) c = O(1) lg = O(5) ll = O(3) qg = O(0)  O(1) ql = O(3)

Rel = O(4) Prg = O(0) Prl = O(0) Scg = O(1) Scl = O(3)

volume fractions prior to achieving steady state. Both issues are discussed below. 2.4.1. Avoiding singularities Oliveira and Issa (2003) pointed out that if the volume fractions of one of the phases goes to zero, the equations become singular. This is especially a problem for the momentum, enthalpy and species equations. They discarded the enthalpy equation and considered a single component flow only and mentioned that prior to reaching this limit the resulting velocities might show large fluctuations. It is inferred that the same sort of fluctuations can be expected for the temperature and mass fractions. Oliveira and Issa (2003) also devised a way to get around this problem. They start by rewriting the equation in non-conservative form. If an arbitrary equation for a variable / is considered, the rewritten equation takes the form

/k

  @ k qk @ k qk uk D/k þ k qk þ ¼ S/ k : @t @z Dt

ð34Þ

The part between brackets is the same as the continuity equation. Oliveira and Issa (2003) considered a two-phase flow without interfacial mass transfer and consequently this term is zero. In the current case the part between brackets should be replaced by the mass source term. Now both sides can be divided by the volume fraction k resulting in

qk

@p þ ui;k Ci;k þ M i;k ; @z

¼ k

15

D/k S/ k /k ¼  C : Dt k k i;k

ð35Þ

Note that both source terms are a function of the volume fraction and would go to zero when the volume fraction goes to zero. Consequently, both terms tend to zero for k ? 0. In that case the equation is no longer singular and the variable reaches its termination value. This approach works only for a fluid that is vanishing and cannot be used for a fluid that is created. In that case / and k are initially zero, which results in a singular equation. Because the liquid volume fraction for the current application is expected to go to zero, this approach will be applied to all equations except the continuity equation. 2.4.2. Bounding of mass and volume fractions The mass and volume fractions are bounded; 0 6 k 6 1 and 0 6 Ya 6 1. Furthermore, they are both subject to a constraint equation given by

X /n ¼ 1;

ð36Þ

n

where / is represents a general variable and n the number of components. For the volume fractions and liquid mass fractions n = 2 and for the gas mass fractions n = 4. Patankar (1980) showed that as long as the coefficients in the linearised algebraic equations are positive, the variable that is solved for is always positive. To ensure that /n is bounded by one, Carver (1982) devised an underrelaxation method where the underrelaxation factor becomes smaller when / approaches one. The variable / is solved for only one of the components, while the other one is found by applying the constraint equation, Eq. (36). However, in the method by Carver (1982) the underrelaxation factor is always larger than zero and consequently does not necessarily guarantee an upper bound of one. Furthermore, although this approach might work for situations where n = 2, boundedness is not guaranteed for n > 2, as is the case for the gas phase. A solution to this was presented by Oliveira and Issa (2003). In their method / is solved for every component. To guarantee an upper bound they normalise / as

16

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

/ /n ¼ P n  ; n /n

ð37Þ

where ⁄ indicates the initial solution of /. When overall converP gence is reached n /n ¼ 1 and thus /n ¼ /n . 2.5. Final set of equations Taking into account the above results in a set of conservation equations that is much simpler to solve and is more stable. The momentum conservation equation for phase k is given by

qk

Duk @p ðui;k  uk ÞCi;k þ Mi;k ¼ þ ; @z Dt k

ð38Þ

and the enthalpy conservation equation by

ap 

Dh Dp ðh  hk ÞCi;k þ Hi;k qk k ¼ þ i;k : Dt Dt k

ð39Þ

The main difference with Eqs. (2) and (3) is caused by ignoring the stress terms. The species conservation is given by

qk

where the velocity correction of the second term on the RHS of Eq. (39) has been ignored. As long as the simulations converge this is justified as at convergence the corrections are zero by definition. Substituting Eq. (42) together with Eqs. (43) and (44) into Eq. (41) results in the pressure-correction equation. After the pressure-correction value has been calculated, the velocities and gas density can be determined from Eq. (42). To ensure stability of the system source terms are split into positive and negative valued groups. Positive source terms are treated explicitly while negative source terms for positive variables are treated implicitly (Patankar, 1980). For an arbitrary variable / the resulting linearised algebraic equations are generally written as

DY k;a Sk;a Y k;a ¼  C : Dt k k i;k

ð40Þ

which is solved for all species and then normalised according to Eq. (37). The continuity equation has not been simplified any further and remains the same as Eq. (1). 3. Solution algorithm In order to solve the system of equations the finite volume method is employed (Ferziger and Peric´, 1999). For this purpose the equations are discretised on a uniform equisized grid. To avoid unwanted decoupling of the pressure and velocity fields a staggered grid arrangement is used, where the scalars are stored at the cell center and the velocities at the cell faces. Cell face fluxes are based on the first order upwind scheme and first-order implicit time stepping is used. As a one-dimensional model is under consideration a tri-diagonal matrix algorithm can be used to solve for each variable. The pressure–velocity coupling is established with a two-phase version of the SIMPLE-algorithm, which is based on the global mass conservation. In this algorithm the gas phase is treated as a compressible phase while the liquid phase is considered to be incompressible although the density can change over time due to a change in temperature and composition. Global mass conservation can thus be written as

 X@ k q @ k qk uk k ¼ 0: þ @t @z k

ð41Þ

! S/ þ v /p ¼ aw /w þ S/ ; /pre p

where S represents the negative part of the source term and S+ the positive part. /prev refers to the value of / from the previous iteration. At the inlet the velocity, temperature, mass fractions and volume fraction for each phase are specified. All values should be non-zero, but can be set to very small values. This makes it possible to perform a simulation where the volume fraction for the gas is so small that the inlet gas mass, momentum and energy flux is negligible. The pressure is specified on the outlet. The solution procedure of the proposed model for each time step is: 1. Solve for the velocities with the momentum conservation Eq. (38) for each phase with a guessed pressure field. 2. Solve the pressure-correction equation based on global mass conservation. 3. Update the velocities, pressure and gas density. 4. Solve for the volume fractions by solving the continuity Eq. (1) for each phase. 5. Solve the species conservation Eq. (40) for each phase. 6. Solve the enthalpy conservation Eq. (39) for each phase. 7. Repeat the steps until convergence has been reached. 4. Validation To validate the model an instrumented catalyst bed has been designed and built. The basic design is similar to the instrumented catalyst bed described in detail by Palmer et al. (2011b), except that the mass flux is higher and there is a greater density of instrumentation on this new bed. A picture of the catalyst bed, including

Initially, the values for qk and uk are calculated from a guessed pressure field p⁄ resulting in guessed values for the density and velocity, qk and uk . A correction, indicated by 0 , needs to be added to the guessed values to obtain the true values

qg ¼ qg þ q0g ; uk ¼ uk þ u0k ; p ¼ p þ p0 :

ð42Þ

Note that q ¼ 0 as the liquid is incompressible. q can be expressed in terms of p0 by means of the equations of state 0 l

q0g ¼

p0 : RT g

0 g

ð43Þ

u0k can be expressed in terms of p0 as well by subtracting the momentum equation, Eq. (27), with velocity uk from the momentum equation with velocity field uk. This results in

qk

Du0k @p0 ¼ ; Dt @z

ð45Þ

ð44Þ Fig. 1. Instrumented catalyst bed used for validation.

17

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

the solenoid valve and standpipes, is shown in Fig. 1. The peroxide flow is from left to right. The catalyst bed has an internal diameter of 16 mm and a length of 96 mm and is made of AISI 316 stainless steel. The inlet is provided by an injector plate with 4 holes with a diameter each of 0.4 mm. The injector plate is connected to a solenoid valve which controls the opening and closing of hydrogen peroxide supply to the catalyst bed. To protect the solenoid valve against damage from thermal soak-back a mass of stainless steel is placed between the injector plate and the solenoid valve which acts as a thermal sink. Upstream of the solenoid valve a mass flow meter is placed (not shown in the picture) and is connected to a propellant tank which is pressurised with nitrogen. The system is designed for a feed pressure up to 20 bar. A retainer plate is located at the downstream side of the bed to keep the catalyst pellets in the catalyst bed. A nozzle with a convergent part only is used to pressurise the bed exit. The throat diameter of the nozzle is 3.55 mm. Three standpipes are located at every 16 mm along the catalyst, which provide access for thermocouples and pressure transducers. At every location the pressure, axial and wall temperature is measured, however, in this study only the pressure and axial temperature are of interest. A heating strip to heat up the catalyst bed to a predetermined temperature is being used prior to firing. An additional set of two standpipes is located between the solenoid valve and injector plate to measure the inlet pressure and temperature. Two other standpipes are placed between the retainer plate and the nozzle to measure the nozzle plenum pressure and temperature. All the transducers are connected to a data acquisition system which measures and records the pressure, temperature and mass flow rate and controls the operating of the solenoid valve and heating strip. For validation the catalyst bed was filled with platinum coated a-alumina cylindrical pellets which had a mean diameter of 3.22 mm and a length of 3.35 mm. The resulting bed to pellet diameter ratio is about 4.3, where the pellet diameter is calculated as the diameter of a sphere with the same volume as the cylindrical pellet. This diameter ratio is so low that the discrete effects of the pellets most likely influence the hydrodynamics. The tuning of the equilibrium constant Kr introduced in Section 2.2.1 and the equivalent radius in the evaporation model introduced in Section 2.2.2 will account for this. This also means that significantly changing the inlet mass flow rate and pellet sizes will require retuning of the model. At the same time from the work of Bey and Eigenberger (1997) can be deduced that the influence of cylindrically shaped pellets on the average void fraction is significantly less than the influence of spherical pellets. During the tests 87.5 wt.% hydrogen peroxide was used, fed by a propellant tank at a fixed pressure of about 17 bar. The catalyst bed was wrapped in an insulation blanket to reduce the heat transfer through the wall as much as possible. Note that in Fig. 1 the catalyst bed is shown without insulation. The measurement frequency was set to 200 Hz and to 1 Hz for temperature measurements. The bed was preheated to 150 °C to reduce the time required to reach steady state. The tank was filled with enough peroxide to run the catalyst bed for about 30 s at which point it was assumed steady state conditions will have been reached. The model described above simulates only what is happening between the injector and retainer plate. The measured values just downstream of the injector and just upstream of the retainer are therefore used as boundary conditions for the numerical model. The outlet pressure boundary condition was set equal to the measured pressure just before the retainer plate, equal to 10.8  105 Pa. At the inlet the gas volume fraction was set to 0.05 and the gas mass fractions to 0.001 for hydrogen peroxide, steam and oxygen and 0.997 for air. The inlet temperature was set equal to the measured peroxide temperature just before filling the pro-

pellant tank. The inlet gas velocity was set to 0.07 m s1 and the liquid velocity was derived from the measured mass flow rate at steady state of 12.6  103 kg s1, which is equivalent to a mass flux of about 62.5 kg m2 s1. By setting these inlet conditions for the gas phase it is assured that a two-fluid flow condition exists at the inlet, but with negligible contribution to the overall inlet mass, momentum and heat flux. Values for the Arrhenius parameters are different for each catalyst material and for each material different researchers have found different values. Instead of choosing a specific catalyst material with corresponding Arrhenius parameters, a set of Arrhenius parameters will be chosen that are representative for typical catalyst for hydrogen peroxide. For liquid thermal decomposition the Arrhenius parameters as determined by Takagi and Ishigure (1985) will be used, who reported for the activation energy EA = 71  103 J mol1 and for the pre-exponential factor A0 = 6.3  105 s1. For thermal decomposition in the gas phase the values found by Hoare et al. (1959) (EA = 200  103 J mol 1 and A0 = 1  1015 s1) will be adopted. For catalytic decomposition a large range of values can be found in the literature. Values for activation energy for liquid and gas will assumed to be EA(l) = 52.5  103 J mol1 and EA = 41.8  103 J mol1 respectively. The former is an average between the value found by Albers et al. (2004) (EA(l) = 50.7  103 J mol1) who investigated manganese oxide catalysts and a research quoted in their work who found a value of EA(l) = 54.5  103 J mol1. In the literature the values of the pre-exponential factor that are quoted appear to vary by at least eight orders of magnitude. This large variation is in part due to the different units that are employed. Catalytic decomposition is a surface reaction for which the unit of the pre-exponential factor is m s1 as opposed to thermal decomposition which is a volumetric process with the unit for the pre-exponential factor of s1. For liquid catalytic decomposition the value found by Albers et al. (2004) will be adopted (A0 = 2.6  104 m s1) No value is reported for A0 for vapour catalytic decomposition. However, Hoare et al. (1959) mentioned that their experiments on homogeneous decomposition a change in reaction rate constant was observed that took place at about 420 °C. They explained this as the transition from catalytic dominated decomposition at low temperatures to thermal dominated decomposition at high temperatures. From this information the pre-exponential factor has been determined to be A0 = 1  101 m s1 while taking into account the reactor surface and volume which was used for the experiment. A summary of the Arrhenius parameters used for the simulations is given in Table 3. All other relevant fluid properties can be found in Appendix A. After the experiments had been performed the catalyst bed was carefully emptied and the number of catalyst pellets was counted. In total there were 406 pellets which, for the dimensions of the pellets and catalyst bed given earlier, results in a void fraction of 0.4262. Using Eq. (5) a void fraction of 0.4215 is predicted which shows very good agreement with the experimental result. The effect of the mesh density is shown in Fig. 2, where the pressure as function of distance from the injector is plotted for difTable 3 Summary of the chosen Arrhenius parameters. Parameter

Value

Unit

A0(l) catalytic decomposition A0(l) thermal decomposition A0(g) catalytic decomposition A0(g) thermal decomposition EA(l) catalytic decomposition EA(l) thermal decomposition EA(g) catalytic decomposition EA(g) thermal decomposition

2.6  104 6.3  105 1  101 1  1015 52.5  103 71  103 41.8  103 200  103

m s1 s1 m s1 s1 J mol1 J mol 1 J mol 1 J mol 1

18

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

ferent grid sizes. A non-dimensional pressure is obtained by dividing by the outlet pressure. Grid independence is shown for grids of 1600 cells and more for the current length of the catalyst bed. All the results hereafter were obtained with this number of cells. The figures shown are the result of steady state simulations. However, when steady state is reached the values are actually undulating around a mean value and is particularly noticeable in the gas temperature. This is caused by small liquid temperature fluctuations near the boiling point. Close to the boiling point, the Spalding mass transfer BM, which drives the evaporation rate, is very sensitive to temperature changes. This results in fluctuating evaporation rates and thus in fluctuating amounts of energy that is consumed during evaporation. For this reason an average was taken over 300 iterations to smooth the fluctuations as much a possible. In that case the standard deviation for the liquid temperature is less than 1.5%. The measured temperature as a function of the distance from the injector together with the simulated gas and liquid temperature is shown in Fig. 3. The distance is again shown in L/D-ratios and the temperature is non-dimensionalised by dividing it by the local boiling temperature, which is a function of the local pressure and composition of the liquid. The measured temperature is below the boiling point up to and including L/D = 2 and gradually rising from the injector. At L/D = 3 a large increase in temperature is shown to a boiling ratio of about 1.8 after which it gradually increases to about 2 at L/D = 5. The temperature then decreases to a ratio of about 1.6 at the end of the catalyst bed. The measured temperature is in fact an average temperature over the measurement time. The averaging process not only smooths the temperature fluctuations within a phase, but also the temperature differences between the phases. Where the temperature is well above the boiling point of the liquid, it may be assumed that only gas is present and thus that the measured temperature is in fact the gas temperature. This assumption is based on two arguments; one which is more mathematical and the other which is more physical in nature. From a mathematical point of view if liquid would be present its influence on the average temperature is so low that it can be regarded as not being present. From a physical point of view if the average temperature is well above the boiling point of the liquid the thermal imbalance between the phases, if a small amount of liquid is present, is large. Consequently the interphase heat transfer will be large which will lead to a rapid increase in liquid temperature, or higher evaporation rates, if the liquid is at the boiling point. Based on this and the measurement result shown in Fig. 3 it can be assumed that liquid is present up to a distance of at least L/D = 2

Fig. 2. Required grid density.

Fig. 3. Temperature ratio of experimental data and numerical results.

from the injector and no liquid is present beyond L/D = 3. The decrease in temperature at the end of the catalyst bed is most likely due to non-adiabatic effects. The simulated liquid and gas temperature are almost the same for temperature ratios below one. In this part of the catalyst bed the gas temperature is mainly the result of decomposition in the liquid phase as can be deduced from the reaction constants and the values of the Arrhenius parameters given in Table 3. The simulated liquid temperature does not exceed a boiling ratio of one, whereas the gas temperature rises to about 2 and stays almost constant in the remainder of the catalyst bed. The figure also shows a large difference between measured and simulated temperature for L/D < 3. This is most likely caused by neglecting the thermal mass and heat transport in the catalyst pellets. This is supported by data shown in Fig. 4, which shows the cell Péclet number for thermal diffusion as a function of the L/D-ratio. Here the Pécelt number is defined as Pe = uDx/b where Dx is the length of the control volume and b the thermal diffusivity defined as qcp/k where k is the thermal conductivity. First of all Fig. 4 shows that thermal diffusion in both fluids can safely be neglected as Peg and Pel, Péclet number of the gas and liquid respectively, is well above 10 except Peg for L/D < 0.5. But as will be shown later,

Fig. 4. Cell Péclet numbers.

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

for this part of the catalyst bed the gas volume fraction is low and consequently the effect of thermal diffusion in the gas phase very small. Pesg and Pesl are the Péclet numbers for the gas and liquid respectively when the pellet thermal diffusivity is taken, with the values for q, cp and k for the pellets as given by Table 4. It shows that for L/D < 3 the thermal diffusion inside the pellets is significant and cannot be neglected. For L/D > = 3, which coincides with the part of the catalyst bed where the boiling ratio is larger than unity, transport by convection is far more important and thermal diffusion can be ignored. The measured and simulated pressure as a function of the distance from the injector in L/D-ratios is shown in Fig. 5 where the pressure is divided by the outlet pressure to make it non-dimensional. Experimental data are plotted with a boxplot with the mean of the measurements superimposed on that. The mean and the median of the measurements are almost identical. The pressure at L/D = 1 shows a peak. However, it is believed that the reason for this is a faulty pressure transducer rather than a physical event taking place at that location. Also shown is the inaccuracy margin of the pressure transducer, which was generated with a MonteCarlo simulation of the inaccuracy of the pressure transducer and data acquisition system. The margin is about 0.9 bar wide. All the measurement data fall within this margin. The solid line is the simulation result. Two regions can be distinguished: a region where the pressure hardly changes and a region where the pressure decreases more rapidly with distance. The pressure hardly changes close to the injector and coincides with the part of the catalyst bed where liquid is present as discussed above. The part of the catalyst bed where the pressure does change coincides with the part that was identified above as the part where only gas is present. A small transition region between the two regions exists. In Section 2.2.1 the equilibrium constant Kr was introduced for catalytic decomposition. It was mentioned that a value has to be found experimentally. The value influences the location at which the liquid temperature reaches a boiling ratio of one and the gas temperature its maximum temperature. The higher the value of Kr, the higher the reaction rates will be, see Eqs. (13) and (14), the closer to the injector this location will be. And the closer this location is to the injector, the larger the pressure drop over the catalyst bed will be as a larger portion of the bed contains gas, which is responsible for most of the pressure drop. The simulations shown in Figs. 3 and 5 were obtained with values of Kr = 3 for the liquid and Kr = 400 for the gas. Fig. 6 shows the location at which the boiling ratio is unity for the first time as a function of the equilibrium constant for the liquid phase. Also shown on the vertical axis is the percent deviation from the measured pressure drop over the catalyst bed for a given equilibrium constant. A similar plot for the gas equilibrium constant is shown in Fig. 7. The vertical axis on the left side shows the distance between the point where the boiling ratio is unity for the first time and the point where the gas temperature reaches 95% of the maximum gas temperature for the first time. The vertical axis on the right side shows the percent deviation from the measured pressure drop over the catalyst bed. Note that the results for the gas phase are not as smooth as for the liquid phase. This has

Table 4 Properties for a-alumina pellets (Munro, 1997). Parameter

Value

Unit

Density Heat capacity Thermal conductivity

3984 755 33

kg m3 J kg1 K1 W m1 K1

19

Fig. 5. Measured and simulated pressure data.

Fig. 6. Influence of equilibrium constant of the liquid phase on pressure and ‘reaction’ distance.

Fig. 7. Influence of equilibrium constant of the gas phase on pressure and ‘reaction’ distance.

to do with the small fluctuations in temperature at the boiling point as explained above. The value for Kr,l should be chosen such that the L/D-ratio where the boiling ratio reaches unity for the first time is close to but larger than 2 for a minimum percentage deviation in pressure. A value of Kr,l = 3.3 is chosen for the liquid equilibrium constant, which

20

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

gives a distance of L/D 2.2 and a pressure overestimation of just under 2%. Based on the measured temperature in Fig. 3, the distance between the point where the boiling ratio is unity for the first time and the point where the gas temperature reaches 95% of the maximum gas temperature for the first time is L/D 0.8. An equilibrium constant for the gas phase of Kr,g = 400 gives a distance of slightly less than 0.8 and an overestimation of the pressure of about 0.2%. As can been seen from Fig. 7 the model is not very sensitive to the value of the gas equilibrium constant, but much more sensitive to the liquid equilibrium constant, see Fig. 6. In Section 2.2.2 the multicomponent evaporation model was presented. The volume equivalent radius rV,a was introduced as a scaling factor to account for the presence of other components at the droplet surface. It was argued that due to the different modes of evaporation in the catalyst bed and due to the consumption of peroxide by decomposition that this factor should be determined experimentally. rV,a has an influence on the maximum gas temperature and the pressure as is shown in Figs. 8 and 9. In both plots the vertical axis on the left side shows the percentage deviation of maximum achieved gas temperature of 965 K as determined with the properties from Appendix A. Here the percentage deviation is calculated as

T D% ¼

T ad  T max;sim ; T ad  T inlet

Fig. 9. Influence of rV;H2 O on pressure and gas temperature.

ð46Þ

where Tad is the adiabatic decomposition temperature, Tmax,sim the maximum simulated temperature and Tinlet the inlet temperature. The axis on the right side shows the percentage deviation from the measured pressure drop over the catalyst bed. The simulations shown in Figs. 3 and 5 were obtained with values of rV;H2 O2 ¼ 0:5 for the liquid and r V;H2 O ¼ 3 for the gas. The values for rV,a should be set such that the maximum simulated temperature is as close as possible to the adiabatic temperature, while keeping the deviation in simulated pressure from measured pressure as small as possible. Based on these restrictions for hydrogen peroxide a value of r V;H2 O2 ¼ 0:5 was set and for water a value of r V;H2 O ¼ 2:5. Substituting the experimentally determined value for the equilibrium constants, Kr, and volume equivalent radius, rV, into the model results in a pressure and temperature profile as shown in Figs. 10 and 11. Fig. 10 shows the temperature normalised to the local boiling temperature together with the experimental data as a function of the distance from the injector. Fig. 11 shows the same for pressure data normalised to the bed exit pressure. The temperature below a boiling ratio of 1 currently cannot be predicted accurately as was explained above, but the sudden in-

Fig. 10. Temperature ratio after tuning the model.

Fig. 11. Pressure data after tuning the model.

Fig. 8. Influence of r V ;H2 O2 on pressure and gas temperature.

crease in temperature and the distance over which this increase takes place shows reasonable agreement with experimental data. Despite the fact that the model assumes adiabatic conditions, the location of the maximum temperature is predicted reasonably well. The overall deviation between measured and simulated pressure drop over the catalyst bed is less than 2.5%, where the measured pressure drop was 1.2 bar absolute. The pressure drop

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

21

gradient is slightly underestimated for the part of the catalyst bed where the boiling ratio is less than unity and slightly overestimated for the rest of the catalyst bed. However, the simulated pressure drop is for most part within the accuracy uncertainty of the measurements. 5. Results and discussion A steady state simulation was performed with the same geometry as the instrumented catalyst bed and initial and boundary conditions as shown in Table 5. The initial concentration of hydrogen peroxide was 87.5% and the initial mass fraction for gaseous hydrogen peroxide, water vapour and oxygen were all set to 0.001 with the rest of the gas being air. Fig. 12 shows the results for the gas and liquid temperature and gas volume fraction. The vertical axis on the left side shows the temperature as a fraction of the maximum temperature which can be expressed in a similar way as Eq. (46)

T k;rel ¼

T k;sim  T k;inlet ; T k;max  T k;inlet

ð47Þ

where Tk,sim is the simulation result of the temperature, Tk,inlet the inlet temperature and Tk,max the local boiling temperature when the liquid phase is considered and the adiabatic decomposition temperature when the gas phase is considered. The vertical axis on the right side shows the gas volume fraction. Both temperatures and gas volume fraction are plotted against the distance from the injector expressed in L/D-ratios. Note that due to a different normalisation of the gas temperature, the result presented in Figs. 3 and 10 differs from the result shown in Fig. 12. The liquid temperature gradually increases from zero until 0.2 from the injector up to L/D = 1.5 after which a sharp increase to about 1 is observed over a distance of just L/D = 0.25. It stays at its boiling point for a distance of about L/D = 0.5 and the gradually decreases to a value of about 0.9. The gas volume fraction shows the same initial behaviour as the liquid temperature. It increases gradually to about 0.3 at L/D = 1.5 from the injector and then a very rapid increase to nearly 1 over a distance of just L/D = 0.25. The gas temperature shows an initial gradual increase to about 0.05 at L/D = 1.5. Then a sharp increase to about 0.3 is observed over a distance of about L/D = 0.25 and then a less sharp, but still significant, increase to 0.95 over a distance of L/D = 0.75. From L/D = 2.5 the increase in gas temperature is gradual until it reaches its maximum temperature at about L/D = 3.5. After that the temperature slowly decreases, despite the fact that the system is assumed to be adiabatic. The reason for this is the heat transfer between the gas and the liquid, which causes the gas temperature to slowly decrease while the liquid temperature reaches a constant temperature for L/D > 4. Note that the gas volume fraction is not exactly one and consequently there is still a small amount of liquid left. This is also shown in Fig. 13 in which the liquid volume fraction

Fig. 12. Temperature and gas volume fraction result.

is plotted on a log scale as a function of the distance from the injector. The decrease in liquid volume fraction is slowed down. As there is still evaporation taking place in the liquid phase, which consumes heat, the liquid is not heated up and an equilibrium is reached between the energy transferred from the gas to the liquid and the energy consumed by evaporation. Figs. 12 and 13 also provide a solid indication that the method to stabilise the scalar equations, as proposed by Oliveira and Issa (2003) and further developed in Section 2.4.1, works. It furthermore shows that even in an adiabatic catalyst bed a maximum gas temperature exists, which is very close but not equal to the adiabatic decomposition temperature, in cases where there is a small amount of liquid left. This result cannot be reproduced in a mixture model as the energy of the mixture is considered and not of each fluid separately. Heat transfer between the fluids is only possible in a two-fluid model. The increase in gas temperature beyond L/D  ratio = 1.75 until the maximum as temperature is reached, is caused by decomposition in the gas phase. Proof of this can be seen in Fig. 14, which shows the gas mass fractions as a function of the L/D-ratio. Note that the mass fraction of air is not shown. At L/D = 1.75 the peroxide vapour mass fraction suddenly increases to about 21% as soon as the liquid temperature is at the boiling point. The corresponding mole fraction of peroxide vapour at this point is about 0.15. Satterfield et al. (1959) conducted a series of ignition tests with hydrogen peroxide vapour at above atmospheric conditions. They showed that for peroxide mole fractions of about 0.20 and more peroxide

Table 5 Geometry, initial and boundary conditions. Parameter

Value

Unit

Catalyst bed length Catalyst bed diameter Pellet shape Pellet length Pellet diameter Inlet liquid mass flux Inlet gas velocity Inlet temperature Gas inlet volume fraction Bed exit pressure

9.6 16 Cylindrical 3.35 3.22 50 0.07 293 0.05 10

cm mm – mm mm kg m2 s1 m s1 K – bar

Fig. 13. Liquid volume fraction as a function of the distance from the injector.

22

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

umix ¼

Fig. 14. Gas mass fractions results.

will ignite spontaneously in the vicinity of a suitable catalyst. This limit remained constant at higher pressures. The peroxide vapour mole fraction of 0.15 from the simulation stays well below this ignition limit and explains why a catalyst bed decomposing 87.5 wt.% peroxide does not explode. It should be noted that Satterfield et al. (1959) did not conduct any experiments at a pressure over 5 bar, but that the experiments at 3 bar and higher suggested that the critical mole fraction would stay constant. As soon as the gas volume fraction is nearly 1, the peroxide vapour mass fraction reaches is maximum and subsequently decreases until it is almost zero at L/D  ratio = 2.5. At this point the peroxide is practically fully decomposed. When the gas volume fraction is nearly one, a negligible amount of liquid hydrogen peroxide will evaporate. As decomposition in the gas phase continues, peroxide is consumed and energy is released, which causes the peak in peroxide vapour mass fraction. Also shown in Fig. 14 is that initially the gas phase composition is mainly oxygen indicating that the majority of the mass transfer is caused by decomposition. At the boiling point a lot of steam and peroxide vapour is generated due to evaporation causing a sharp decrease in the oxygen mass fraction and a sharp increase in the hydrogen peroxide vapour and steam mass fraction. In Section 2.2.3 a model was introduced to describe the pressure drop resulting from the interaction between the fluids and the fluids with the catalyst bed. Traditionally, the Ergun-equation is used for this (Bird et al., 2002) such as was done by Pasini et al. (2010) and written as

  1  f qu2 Dp A ¼ þB ; L Re 3f dp

ð48Þ

g qg ug þ l ql ul qmix

:

ð50Þ

qmix and umix were then substituted into Eq. (48). The result of this was then divided by the pressure gradient as was computed by the two-fluid model. The resulting pressure gradient for mixture models for the original Ergun-equation, the Ergun-equation with the modification proposed by Macdonald et al. (1979) and the modification proposed by Tallmadge (1970) is shown in Fig. 15. For the parameter B for the Macdonald modification the value B = 1.8 was used, corresponding to smooth pellets. In this figure the gradients are divided by the gradient as computed by the two-fluid model, a pressure gradient ratio higher than 1 means that a higher pressure drop is achieved. First of all it is evident that original Ergun-equation as well as the modification proposed by Macdonald et al. (1979) result in very much higher pressure gradients than the extension proposed by Tallmadge (1970). Especially when the gas temperature approaches its maximum, the Tallmadge equation gives a pressure gradient which is at least a factor 2 lower. Also shown is that the Macdonald modification results in higher pressure gradients than the original Ergun-equation. Note that in this plot smooth pellets were assumed. For less smooth pellets the value for B would increase resulting in even higher gradients. The second point that becomes clear in this figure is that all mixture models give a higher pressure gradient than the two-fluid model except for a small region around L/D = 1.75. Only in the region where the temperature reaches the boiling point for the first time the ratio, and thus the pressure gradient, is lower than unity. The reason for this difference can be explained with Figs. 13 and 16. The latter shows the gas and liquid phase velocity as computed by the two-fluid model divided by the mixture velocity umix determined by Eq. (50). From the Result 1 is subtracted to enhance the readability: values below zero indicate a velocity less than the mixture velocity and values above zero indicate values higher than the mixture velocity. The relative gas velocity in Fig. 16 is about 0.17, which means that it is about 17% higher than the mixture velocity. The relative liquid velocity is about 0.92% or 92% lower than umix. As is shown in Fig. 13, there is still a small amount of liquid left in the catalyst bed. Consequently, part of the pressure gradient is due to the (higher than umix) gas velocity and part due to the (lower than umix) liquid velocity. About 90% of the total mass transfer between the gas and liquid phase takes place in a very narrow section of the catalyst bed: for the current simulation it is 1.75 < L/D < 2 as is shown in Fig. 17. The solid line in this figure indicates the cumulative mass transfer as a function of the distance to the injector. Only a small amount, less

where in the original Ergun-equation A = 150 and B = 1.75. Macdonald et al. (1979) tested this equation against several sets of experimental data and modified the constants A and B to A = 180 and B = 1.8 4 where the exact value of B is dependent on how smooth the particles are. As was pointed out in Section 2.2.3 applicability of the Ergun-equation is limited to lower Reynolds numbers, typically 101 < Re < 103. From Table 2 can be seen that the expected Re number is of order 104. Tallmadge (1970) extended the equation to included higher Reynolds numbers, see Eq. (25). To investigate the influence of the different pressure models on the overall pressure drop as well as to investigate the difference between mixture and two-fluid models, the mixture pressure gradient has been computed based on the results for the two-fluid model. For this purpose first the mixture density, qmix, and mixture velocity, umix, were determined by

qmix ¼ g ug þ l ul ;

ð49Þ

Fig. 15. Pressure gradient ratio for several mixture pressure drop models.

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

23

total mass transfer and the cumulative total mass transfer. This plot supports the argument given earlier that decomposition is responsible for the mass transfer close to the inlet. Based on these plots three regions can be distinguished inside the catalyst bed. These are:

pre-boiling region: from the injector to the point where the liquid reaches the boiling point for the first time,

rapid conversion region: from the point where the liquid reaches the boiling point for the first time to the point where the gas volume fraction is more than 0.99,

dry-out region: from the point where the gas volume fraction is more than 0.99 to the outlet.

Fig. 16. Gas and liquid velocity relative to the mixture velocity.

In the pre-boiling region decomposition accounts for most of the mass transfer between the gas and liquid, but accounts for only a small fraction of the total mass transfer in the catalyst bed. Despite that, the gas volume fraction increases significantly. In the rapid conversion region about 90% of the total mass transfer takes place mainly by means of evaporation. This is also the region where there is the largest liquid temperature increase and where the gas volume fraction becomes nearly unity. The dry-out region is characterised by a significant increase in gas temperature and only a small amount of mass transfer. The three regions identified above can also be identified in the works by Zhou and Hitt (2003) and Pasini et al. (2010). The results presented by Zhou and Hitt (2003) also show a very rapid increase in temperature from well below the boiling point to a temperature far above it. However, any information on the gas and liquid volume fractions cannot be obtained from their model. Also the relative contribution of each transfer mode cannot be obtained. The results presented by Pasini et al. (2010) show a gradual increase in temperature and an almost linear variation in pressure drop. Similar to the work of Zhou and Hitt, no information is available on volume fractions. This shows once more the advantages of the two-fluid model developed in this paper.

Fig. 17. Cumulative relative mass transfer.

6. Conclusions

Fig. 18. Relative contribution of decomposition to mass transfer.

than 5%, takes place between the inlet and L/D = 1.75. Also shown is the relative contribution of decomposition, mass transfer due to formation of oxygen during decomposition of the liquid phase, and evaporation. Evaporation accounts for the majority of the mass transfer: about 90%. It should be noted that despite it being the major mode of mass transfer, it becomes important only after the liquid approaches the boiling point. This is shown in Fig. 18, in which the ratio is shown of the cumulative decomposition to the

In this paper a two-fluid multicomponent flow model for packed beds has been described and validated. A dimensional analysis has been performed to justify the simplifications of the model. The calculation procedure for the mass and volume fractions is such that a lower bound of zero and an upper bound of unity is guaranteed. The model was applied to a decomposing flow of highly concentrated hydrogen peroxide in a pellet-based catalyst bed. For this purpose models have been presented describing the mass, momentum and heat transfer of both multicomponent phases. An instrumented catalyst bed was used for tuning of the model and validation against experimental results. Overall a reasonable agreement between simulations and experimental data was found. It is shown that the lack of a thermal description of the catalyst material in the current model causes a large difference between the simulated and measured temperature in the upstream region of the catalyst where gas and liquid are both present. However, these differences are not present in the region of the catalyst bed dominated by the gas phase. It was shown that a two-fluid model gives a much more reliable estimation of the pressure drop over the catalyst bed than a mixture model. It was also shown that a two-fluid formulation in general reveals features in the catalyst bed that are not reproducible by mixture models due to the lack of a description of interphase energy transport. Finally, it was shown that the catalyst bed can be subdivided into three regions: pre-boiling region, rapid conversion region and dry-out region. The pre-boiling region is characterised by a rel-

24

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

atively large increase in gas volume fraction driven by decomposition in the liquid phase, but hardly any mass transfer. The rapid conversion region is characterised by a large increase in volume fraction to almost unity and a high liquid to gas mass conversion rate, primarily caused by evaporation. The dry-out phase is characterised by a small amount of mass transfer at a sharply increasing gas temperature. Most of the conversion is achieved by evaporation. These three regions can also be identified in the works of other researchers.

Table A.8 Constants for the density of water. Constant

Value

A B C D E F G

0.9998396 18.224944  103 7.92221  106 55.44846  109 149.7562  1012 393.2952  1015 18.159725  103

Acknowledgement The research leading to these results has received funding from the European Community Seventh Framework Programme (FP7/ 2007-2013) under Grant agreement no 218819 and was conducted under the GRASP (GReen Advanced Space Propulsion) Project. The authors would also like to thank Mr. Palmer, a research student at the University of Southampton, for the design of the instrumented catalyst bed and assisting with the experiments with which validation has been carried out.

Table A.9 Specific heat at constant pressure for liquid hydrogen peroxide (McCormick, 1965) and air. Component

cp (J mol1 K1)

H2O2(l) Air

89.377 29.086

A.4. Specific heat and enthalpies The specific heats in J mol1 K1 for hydrogen peroxide vapour, water, steam and oxygen are given by a polynomial of the form

Appendix A. Fluid properties A.1. Molar masses

cp ¼ A þ B  t þ C  t 2 þ D  t 3 þ See Table A.6.

E ; t2

ðA:3Þ

where t = Ta/1000. The values for the constants are given in Table A.11 below. For liquid hydrogen peroxide and air the specific heat is given by a constant, see Table A.9. The specific enthalpy can be found by integrating the specific heat over the temperature. Applying this to Eq. (A.3) and using a reference temperature of 298.15 K gives

A.2. Enthalpy of formation See Table A.7. A.3. Liquid density

ha ¼ The density of pure liquid hydrogen peroxide and water is temperature dependent. For hydrogen peroxide it is given by (Tsykalo and Tabachnikov, 1966)

298:15

þ G  T l;C Þ  1000:

ðA:2Þ

where Tl,C is the liquid temperature in degrees Celsius and the coefficients are given in Table A.8.

Table A.6 Molar masses (NIST, 2005).

where ha is given in J mol1.

DHfg ¼ DHfg;T b

Table A.7 Enthalpy of formation (NIST, 2005).

1

)



1  Tr 1  T br

q ðA:5Þ

:

In this equation DHfg;T b is the latent heat at the normal boiling point as given in Table A.10, Tr the reduced temperature and Tbr the reduced normal boiling point. The reduced temperature is defined as:

Tr ¼ 34.0147 18.0153 31.9988 28.97

ðA:4Þ

The enthalpy of evaporation for hydrogen peroxide and water is a function of temperature and can be approximated by the Watson equation (Li et al., 1997):

ql;H2 O ¼ ðA þ B  T l;C þ C  T 2l;C þ D  T 3l;C þ E  T 4l;C þ F  T 5l;C Þ=ð1

H2O2 H2O O2 Air

1 1 1 E cp dT ¼ A  t þ B  t2 þ C  t3 þ D  t4  þ F; 2 3 4 t

ðA:1Þ

and for water by (Perry, 1950)

Molar mass (g mol

Ta

A.5. Enthalpy of evaporation

ql;H2 O2 ¼ 1597 þ 0:0784T l  0:00197T 2l ;

Component

Z

T : Tc

ðA:6Þ

Tc is the critical temperature and, together with the normal boiling points, is further specified in Table A.10 (Manatt and Manatt, 2004). The parameter q can be calculated according to Fishtine by (Li et al., 1997)

Table A.10 Latent heat properties (Perry, 1950).

Component

DHf J mol1

Species

DHfg (J mol1)

Tb (K)

Tc (K)

H2O2 H2O

187.86  103 285.83  103

H2O2 H2O

42969.68 40706.136

426.305 373.15

739.5 647.3

25

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

8 for T br < 0:57 > < 0:30 q ¼ 0:740T br  0:116 for 0:57 6 T br 6 0:71 > : 0:41 for T br > 0:71:

ðA:7Þ

A.6. Vapour pressure The vapour pressure in mmHg of pure hydrogen peroxide and pure water as a function of the temperature of the liquid is generally given by the following relation

B ¼ A þ þ C  log10 T l þ D  T l þ E  T 2l þ F  T 3l þ G Tl

log10 pv ap

 T 4l :

ðA:8Þ

The values for the constants are given in Table A.12.

Bi ¼ C 0i þ

C 1i : 1 þ eC 2i ðTC 3i Þ

ðA:12Þ

B0 consists of 4 different regions. For temperatures between 273.150 and 317.636 K the curve is described by a Lorentzian function of the same form as given in Eq. (A.11). For temperatures between 348.222 and 391.463 K the curve takes the form of a second order polynomial:

B0 ¼ P 0 þ P1 T l þ P2 T 2l :

ðA:13Þ

For temperatures between 317.636 and 348.222 K B0 is the average of the Lorentzian curve given by Eq. (A.11) and the second order polynomial given by Eq. (A.13). And finally, for temperature above 391.463 up to 433.150 K B0 is a constant with a value of B0 =  612.9613. All the coefficients mentioned in Eqs. A.11,A.12 and A.13 are given in Table A.13 below together with the temperature range over which they are applicable.

A.7. Activity coefficients A.8. Binary diffusion coefficient The activity coefficient for hydrogen peroxide and water is expressed in terms of the liquid mole fraction and Redlich–Kister parameters (Radl et al., 2009). For hydrogen peroxide it is

ln cH2 O2 ¼

x2H2 O RT l

 ½B0 ðT l Þ þ B1 ðT l Þ  ð3  4xH2 O Þ þ B2 ðT l Þ

ð1  2xH2 O Þð5  6xH2 O Þ þ B3 ðT l Þ  ð1  2xH2 O Þ2 ð7  8xH2 O Þ;

The binary diffusion coefficients for the gas phase species are assumed to be similar and described sufficiently accurate by the binary diffusion coefficient of water vapour in air provided by Bolz and Tuve (1973) (see Table A.11,A.12,A.13).

Da;g ¼ 2:775  106 þ 4:479  108 T g þ 1:656  1010 T 2g :

ðA:9Þ

ðA:14Þ

and for water

ln cH2 O ¼

1  x2H2 O RT l

Appendix B. Interfacial area

 ½B0 ðT l Þ þ B1 ðT l Þ  ð1  4xH2 O Þ þ B2 ðT l Þ

ð1  2xH2 O Þð1  6xH2 O Þ þ B3 ðT l Þ  ð1  2xH2 O Þ2 ð1  8xH2 O Þ: ðA:10Þ The parameters Bi are temperature dependent. B1 is given by a Lorentzian function of the form

B1 ¼ C 0 þ 

p

C1C2 C 22

þ ðT l  C 3 Þ

2

:

ðA:11Þ

The parameters B2 and B3 are described by sigmoid function of the form

The total surface area of the catalyst material is obtained by modelling the catalyst material as a long slender rod with a diameter equal to the diameter of the catalyst particle, dp. The length of this rod, Lrod, is equal to the total volume of the catalyst material divided by the cross sectional area of the pellet

Lrod ¼

s V bed : pð0:5dp Þ2

ðB:1Þ

The part of the rod that is covered with liquid, Lwet, is proportional to the ratio of the liquid volume fraction to the total void space in the catalyst bed. This can be expressed as

Table A.11 Specific heat constants (NIST, 2005). Constant

H2O2(g)

H2O(l)

H2O(g)

O2(T 6 700 K)

O2(T > 700 K)

A B C D E F

34.25667 55.18445 35.15443 9.087440 0.422157 13.8034

203.606 1523.29 3196.413 2474.455 3.855326 29.2826

30.09200 6.832514 6.7934535 2.534480 0.082139 9.0546

31.32234 20.23531 57.86644 36.50624 0.007374 8.903471

30.03235 8.772972 3.988133 0.788313 0.741599 11.32468

Table A.12 Vapour pressure constants (Manatt and Manatt, 2004). Constant

H2O2 (T < 363.15 K)

H2O2 (T P 363.15 K)

H2O

A B C D E F G

24.8436 3511.54 4.61453 3.60245  103 7.73423  106 1.78355  108 2.27008  1013

38.8572 3627.72 11.2133 4.74132  103 0 0 0

19.389127 2861.9133 3.2418662 1.0799994  104 7.9189289  106 1.5411774  108 8.1926991  1012

26

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28 Table A.13 Coefficients for the Redlich–Kister parameters (Manatt and Manatt, 2004). Parameter

Temperature range

Curve type

Constants

B0

273.150–317.636 K

Lorentzian

317.636–348.222 K

Average of above Lorentzian and 2nd order polynomial

C0 =  999.8830 C1 = 2499.584 C2 = 8.261924 C3 = 327.4487 P0 = 17418.34 P1 = 109.9125 P2 = 0.1663847

348.222–391.463 K

2nd order polynomial

P0 = 6110.401 P1 = 28.08669 P2 = 0.03587408 B0 = 612.9613

391.463–433.150 K

Constant

B1

273.150–433.150 K

Lorentzian

C0 = 126.7385 C1 = 2558.776 C2 = 12.33364 C3 = 343.1050

B2

273.150–433.150 K

Sigmoid

C02 = 63.18354 C12 = 149.9278 C22 = 0.4745954 C32 = 348.1642

B3

273.150–433.150 K

Sigmoid

C03 = 59.42228 C13 = 199.2644 C23 = 0.8321514 C33 = 346.2121

Appendix C. Evaporation model The total evaporation rate is given by (Brenn et al., 2007)

_ ev ap ¼ 2p m

X rV;a qg Da;g Sha lnð1 þ BM;a Þ;

ðC:1Þ

a

where the subscript g refers to the gas phase. Note that a negative value for the evaporation rate denotes evaporation. BM;a is the Spalding mass transfer number for component a and is defined as Fig. B.19. Schematic representation of the liquid distribution around the rod.

Lwet ¼ l Lrod :

ðB:2Þ

The interfacial area of the liquid and the solid, Asl, is now equal to the surface area of a cylinder with a diameter of dp and a length of Lwet

Asl ¼ 4l

s V bed dp

ðB:3Þ

:

The remainder of the rod is in contact with the gas phase of which the interfacial area, Asg, can be obtained in a similar way

Asg ¼ 4g

s V bed dp

sffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 1 :

s

s

ðC:2Þ

Sh⁄ is the modified Sherwood number and takes into account the effect of a high mass transfer rate on the evaporation process. It is defined as

Sh ¼ 2 þ

Sh0  2 : FðBM Þ

ðC:3Þ

The Sherwood number Sh0 is calculated as a function of the Reynolds and Schmidt number

Sh0 ¼ 2 þ 0:552Re1=2 Sc1=3 :

ðC:4Þ

ðB:4Þ

ðB:5Þ

The surface area, Alg, can now be worked out and in combination with Eq. (B.3) it can be written as

sffiffiffiffi 1 : Alg ¼ Asl

Y a;s  Y a;1 ; 1  Y a;s

The function F(B) in Eq. (C.3) is generally given by

:

Evaporation is assumed to take place at the gas–liquid surface only, as indicated in Fig. B.19. To determine the surface area it is assumed that the liquid around the wetted part of the rod is distributed homogeneously. The thickness, t, of the liquid layer is such that the volume of the liquid layer is equal to the total volume of liquid in the catalyst bed. After some rearrangement this can be written as

t ¼ 0:5dp

BM;a ¼

ðB:6Þ

FðBÞ ¼ ð1 þ BÞ0:7

lnð1 þ BÞ : B

ðC:5Þ

The heat transfer rate of the droplet is given by

h_ l ¼

  N X _ A cp;A ðT 1  T s Þ m  DHfg ðT s Þ ; BTA md A¼1

ðC:6Þ

where cp is the specific heat at constant pressure, md the mass of the droplet, T1 the temperature of the gas far away from the droplet surface and Ts the temperature of the droplet surface. BTA is given by

BTA ¼ ð1 þ BMA ÞwA  1;

ðC:7Þ

where wA is calculated by

wA ¼

cp;A ShA 1 : cpg NuA Le

ðC:8Þ

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

Le is the Lewis number calculated as Le = Sc/Pr and NuA is the modified Nusselt number, which is calculated in a similar way as the modified Sherwood number

NuA ¼ 2 þ

Nu  2 : FðBTA Þ

ðC:9Þ

The Nusselt number Nu is calculated in the same way as the Sherwood number, see Eq. (C.4), with the Schmidt number substituted by the Prandtl number. Note that an iteration loop between Eqs. (C.8) and (C.9) is required. Appendix D. Pressure drop model The model by Sorokin (2008) considers the pressure drop by each phase separately. In the current model the pressure drop for a single pase is given by the Tallmadge equation (Tallmadge, 1970)

  Dp 5=6 11=6 1=6 ¼ K ll ul þ gql ul ll ; L l

ðD:1Þ

  Dp 11=6 1=6 ¼ K lg ug þ gq5=6 lg ; g ug L g

ðD:2Þ

where

K¼

150 2 dp

2s ð1  s Þ3

and g ¼

4:2 7=6 dp

s7=6 : ð1  s Þ3

The combined pressure drop for a two-phase flow through a packed bed is then calculated by

    Dp Dp ¼ W2 ; L tp L k k

ðD:3Þ

where subscript k refers to either the liquid or gas phase. The parameter W for a liquid and a gas is defined as respectively

W2l ¼ 1 þ

C 1 þ ; X X2

W2g ¼ 1 þ CX þ X 2 :

ðD:4Þ ðD:5Þ

The parameter X is the square root of the ratio of the liquid pressure drop and gas pressure drop in a packed bed and defined as

X¼

s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ffi Dp Dp = : L l L g

ðD:6Þ

The parameter C is defined as

gql r C¼ K ll Jdp

!0:3

;

ðD:7Þ

where r represents the surface tension of the liquid and J is the mass flux. References Abramzon, B., Sirignano, W., 1989. Droplet vaporization model for spray combustion calculations. Int. J. Heat Mass Transfer 32, 1605–1618. Albers, R., Houterman, M., Vergunst, T., Grolman, E., Moulijn, J., 2004. Novel monolithic stirred reactor. AIChE J. 44, 2459–2464. Atkins, P., De Paula, J., 2006. Atkins Physical Chemistry, uk ed. Oxford University Press. Bellan, J., Summerfield, M., 1978. Theoretical examination of assumptions commonly used for the gas phase surrounding a burning droplet. Combust. Flame 33, 107–122. Bey, O., Eigenberger, G., 1997. Fluid flow through catalyst filled tubes. Chem. Eng. Sci. 52, 1365–1376. Bird, R., Stewart, W., Lightfoot, E., 2002. Transport Phenomena, second ed. John Wiley & Sons, USA. Bolz, R., Tuve, G., 1973. Crc Handbook of Tables for Applied Engineering Science (edn).

27

Bonifacio, S., Russo Sorge, A., 2006. Modelling hydrogen peroxide decomposition in monolithic beds. In: Proceedings of 3rd International Conference on Green Propellant for Space Propulsion and 9th International Hydrogen Peroxide Propulsion Conference. SP-635. ESA. Brenn, G., Deviprasath, L., Durst, F., Fink, C., 2007. Evaporation of acoustically levitated multi-component liquid droplets. Int. J. Heat Mass Transfer 50, 5073– 5086. Carver, M., 1982. A method of limiting intermediate values of volume fraction in iterative two-fluid computations. J. Mech. Eng. Sci. 24, 221–224. Chou, S., Huang, C., 1999. Decomposition of hydrogen peroxide in a catalytic fluidized-bed reactor. Appl. Catal. A: Gen. 185, 237–245. Corpening, J., Heister, S., Anderson, W., Austin, B., 2006. Thermal decomposition of hydrogen peroxide, Part 2: modeling studies. J. Propul. Power 22, 996. Fanelli, E., Turco, M., Russo, A., Bagnasco, G., Marchese, S., Pernice, P., Aronne, A., 2011. Mnox/zro2 gel-derived materials for hydrogen peroxide decomposition. J. Sol–Gel Sci. Technol. 60, 426–436. Ferziger, J., Peric´, M., 1999. Computational Methods for Fluid Dynamics. Springer Verlag. Hatsopoulos, G., Keenan, J., 1981. Principles of General Thermodynamics. Krieger. Helms, W., Mok, J., Sisco, J., Anderson, W., 2002. Decomposition and vaporization studies of hydrogen peroxide. In: 38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Indianapolis, USA. No. AIAA 2002-4028. Hoare, D., Protheroe, J., Walsh, A., 1959. The thermal decomposition of hydrogen peroxide vapour. Trans. Faraday Soc. 55, 548–557. Ishii, M., Hibiki, T., 2006. Thermo-Fluid Dynamics of Two-Phase Flow. Springer. Johnson, C., Anderson, W., Ross, R., Lyles, G., 2000. Catalyst bed instability within the usfe h2o2/jp-8 rocket engine. In: 36th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Huntsville, Alabama. No. AIAA 2000-3301. Kappenstein, C., 2008. Contribution to the Green Propulsion Workshop. Crete. Li, P., Liang, Y., Ma, P., Zhu, C., 1997. Estimations of enthalpies of vaporization of pure compounds at different temperatures by a corresponding-states groupcontribution method. Fluid Phase Equilibr. 137, 63–74. Macdonald, I., El-Sayed, M., Mow, K., Dullien, F., 1979. Flow through porous mediathe Ergun equation revisited. Ind. Eng. Chem. Fund. 18, 199–208. Manatt, S., Manatt, M., 2004. On the analyses of mixture vapor pressure data: the hydrogen peroxide/water system and its excess thermodynamic functions. Chem. – A Eur. J. 10, 6540–6557. McCormick, J., 1965. Hydrogen Peroxide Rocket Manual. FMC Corporation. Miller, R., Harstad, K., Bellan, J., 1998. Evaluation of equilibrium and nonequilibrium evaporation models for many-droplet gas–liquid flow simulations. Int. J. Multiphase Flow 24, 1025–1055. Mol, J.C., Ertl, G., Knötzinger, H., Schüth, F., Weitkamp, J., 2008. Handbook of Heterogeneous Catalysis 5, 2647–2676. Munro, M., 1997. Evaluated material properties for a sintered alpha-alumina. J. Am. Ceram. Soc. 80, 1919–1928. Musker, A., Roberts, G., Chandler, P., Grayson, J., Holdsworth, J., 2004. Optimisation study of a homogeneously-catalysed HTP rocket engine. In: ESA Special Publication. vol. 557. pp. 9. NIST, June 2005. Nist Chemistry Webbook. . Oehmichen, T., Datsevich, L., Jess, A., 2010. Influence of bubble evolution on the effective kinetics of heterogeneously catalyzed gas/liquid reactions. Part I: reactions with gaseous products. Chem. Eng. Technol. 33, 911–920. Oliveira, P., Issa, R., 2003. Numerical aspects of an algorithm for the Eulerian simulation of two-phase flows. Int. J. Numer. Methods Fluids 43, 1177–1198. Palmer, M., Musker, A., Roberts, G., 2011a. Design, build and test of a 20N hydrogen peroxide monopropellant thruster. In: 47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, San Diego, California. No. AIAA 20115697. Palmer, M., Musker, A., Roberts, G., 2011b. Experimental assessment of heterogeneous catalysts for the decomposition of hydrogen peroxide. In: 47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, San Diego, California. No. AIAA 2011-5695. Pasini, A., Torre, L., Romeo, L., Cervone, A., D’Agostino, L., 2008. Testing and characterization of a hydrogen peroxide monopropellant thruster. J. Propul. Power 24, 507–515. Pasini, A., Torre, L., Romeo, L., Cervone, A., D’Agostino, L., 2010. Reduced-order model for H2O2 catalytic reactor performance analysis. J. Propul. Power 26, 446– 453. Patankar, S., 1980. Numerical Heat Transfer and Fluid Flow. Taylor & Francis. Pearson, N., Pourpoint, T., Anderson, W., 2003. Vaporization and decomposition of hydrogen peroxide drops. In: 39th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Huntsville, Alabama. No. AIAA 2003-4642. Perry, J., 1950. Chemical engineers’ handbook. J. Chem. Ed. 27, 533. Radl, S., Ortner, S., Sungkorn, R., Khinast, J., 2009. The engineering of hydrogen peroxide decontamination systems. J. Pharm. Innov. 4, 51–62. Sorokin, V., 2008. Calculation of two-phase adiabatic flow in a pebble bed. High Temp. 46, 432–434. Satterfield, C., Feakes, F., Sekler, N., 1959. Ignition limits of hydrogen peroxide vapor at pressures above atmospheric. J. Chem. Eng. Data 4, 131–133. Takagi, J., Ishigure, K., 1985. Thermal decomposition of hydrogen peroxide and its effect on reactor water monitoring of boiling water reactors. Nucl. Sci. Eng.; (United States) 89. Tallmadge, J., 1970. Packed bed pressure drop an extension to higher reynolds numbers. AIChE J. 16, 1092–1093. Tsykalo, A., Tabachnikov, A., 1966. Density, viscosity, and bond energy of molecules in aqueous hydrogen peroxide solutions. Theor. Exp. Chem. 2, 604–605.

28

R.-J. Koopmans et al. / International Journal of Multiphase Flow 57 (2013) 10–28

Wernimont, E., Ventura, M., Garboden, G., Mullens, P., 1999. Past and present uses of rocket grade hydrogen peroxide. In: International Hydrogen Peroxide Propulsion Conference, West Lafayette, USA. Zeigarnik, Y., Kalmykov, I., 1985. Experimental study of the flow resistance of porous structures in the adiabatic motion of steam-water mixtures. High Temp.(Engl. Transl.); (United States) 23.

Zhou, X., Hitt, D., 2003. One-dimensional modeling of catalyzed H2O2 decomposition in microchannel flows. In: 33rd AIAA Fluid Dynamics Conference and Exhibit. No. AIAA 2003-3584. Zhou, X., Hitt, D., 2005. Numerical Modeling of Monopropellant Decomposition in a Micro-Catalyst Bed (AIAA 2005-5033).