The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments

G Model

ARTICLE IN PRESS

NED-7927; No. of Pages 15

Nuclear Engineering and Design xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments Pratap Sathiah c,∗ , Ed Komen a , Dirk Roekaerts b a

Nuclear Research and Consultancy Group – NRG, P.O. Box 25, 1755 ZG Petten, The Netherlands Delft University of Technology, P.O. Box 5, 2600 AA Delft, The Netherlands c Shell Global Solutions Ltd., Brabazon House, Concord Business Park, Threapwood Road, Manchester M220RR, United Kingdom b

h i g h l i g h t s • • • •

A CFD based method proposed in the previous article is used for the simulation of the effect of CO2 –He dilution on hydrogen deflagration. A theoretical study is presented to verify whether CO2 –He diluent can be used as a replacement for H2 O as diluent. CFD model used for the validation work is described. TFC combustion model results are in good agreement with large-scale homogeneous hydrogen–air–CO2 –He experiments.

a r t i c l e

i n f o

Article history: Received 29 October 2013 Received in revised form 28 May 2014 Accepted 30 May 2014

a b s t r a c t Large quantities of hydrogen can be generated and released into the containment during a severe accident in a PWR. The generated hydrogen, when mixed with air, can lead to hydrogen combustion. The dynamic pressure loads resulting from hydrogen combustion can be detrimental to the structural integrity of the reactor safety systems and the reactor containment. Therefore, accurate prediction of these pressure loads is an important safety issue. In our previous article, a CFD based method to determine these pressure loads was presented. This CFD method is based on the application of a turbulent flame speed closure combustion model. The method was validated against three uniform hydrogen–air deflagration experiments with different blockage ratio performed in the ENACCEF facility. It was concluded that the maximum pressures were predicted within 13% accuracy, while the rate of pressure rise dp/dt was predicted within about 30%. The eigen frequencies of the residual pressure wave phenomena were predicted within a few %. In the present article, we perform additional validation of the CFD based method against three uniform hydrogen–air–CO2 –He deflagration experiments with three different concentrations of the CO2 –He diluent. The trends of decrease in the flame velocity, the intermediate peak pressure, the rate of pressure rise dp/dt, and the maximum value of the mean pressure with an increase in the CO2 –He dilution are captured well in the simulations. From the presented validation analyses, it can be concluded that the maximum value of the mean pressures and the intermediate peak pressures were predicted respectively within 12 and 29% accuracy, while the rate of pressure rise dp/dt was typically underpredicted within 15–90%. The eigen frequencies of the residual pressure wave phenomena were predicted within 6%. It was overall concluded that the current model predicts the considered ENACCEF experiments well. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Large quantities of hydrogen can be generated and released into the containment during a severe accident in a Pressurized Water

∗ Corresponding author. Tel.: +447580477607. E-mail address: [email protected] (P. Sathiah).

Reactor (PWR). When mixed with air in the containment, the hydrogen can form a flammable or explosive gas mixture. Upon ignition, a combustion process is initiated that may damage relevant safety systems and that may even cause a threat for the integrity of the containment walls. The most recent examples of hydrogen explosion took place at the Fukushima Daiichi reactors in Japan (Fukushima, 2011). These examples confirmed that the release, distribution, and possible

http://dx.doi.org/10.1016/j.nucengdes.2014.05.042 0029-5493/© 2014 Elsevier B.V. All rights reserved.

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

2

Nomenclature CFD UDF LWR PAR PMT BR AMR AICC CFL CNRS D d dp/dt p t Xsteam B A1 A2 A3 A4 A5 A6 D1 C Tref Pref T P c I ω  b Sl Z L ı

Computational Fluid Dynamics user defined functions Light Water Reactor Passive Autocatalytic Recombiner photomultiplier tubes blockage ratio adaptive mesh refinement Adiabatic Isochoric Complete Combustion Courant–Friedrichs–Lewy Centre National de la Recherche Scientifique diameter of pipe [m] diameter of baffle [m] slope of pressure rise with time [Pa/s] pressure [Pa] time [s] steam concentration model constant model constant model constant model constant model constant model constant model constant model constant model constant reference temperature [K] reference pressure [pa] temperature [K] pressure [pa] speed of sound [m/s] acoustic intensity [N/m s] eigen frequency [Hz] equivalence ratio diluent mole fraction burnt gas density [kg/m3 ] laminar flame speed [m/s] acoustic impedance [kg/m2 s] length of the pipe [m] acoustic fluctuating quantities

combustion of hydrogen are still key safety issues for Light Water Reactors (LWRs). Passive Autocatalytic Recombiners (PARs) are widely used as a mitigation measure to reduce the risk of possible hydrogen combustion. Despite the installation of PARs, it has been generally recognized that the temporary existence of flammable gas clouds cannot be fully excluded during certain postulated accident scenarios (e.g. Bentaib et al., 2010). Therefore, reliable computer modeling is needed to assess the associated residual risk of possible hydrogen deflagrations and to optimally design the hydrogen mitigation systems in order to reduce this risk as far as possible. The computer modeling applied at NRG for hydrogen safety analyses is based on a combined lumped parameter code/Computational Fluid Dynamics (CFD) code modeling approach. The lumped parameter code MELCOR (Summers et al., 1991) or SPECTRA (Stempniewicz, 2009) is used to determine the evolution of the hydrogen distribution within the containment for a relatively large number of postulated accident scenarios. For the most critical scenarios, more refined hydrogen distribution analyses are performed using NRG’s CFD based containment model. For the scenarios for which a fast deflagration cannot be ruled out, hydrogen combustion analyses are performed also with NRG’s CFDbased containment model. The main reason for the need of CFD

modeling for such scenarios can be summarized as follows: The thermal and dynamic pressure loads that occur during a hydrogen deflagration are determined by the turbulent flame acceleration. This turbulent flame acceleration is determined by the amount of turbulence generated locally during the combustion process. In turn, this turbulence generation is highly affected by the obstacles that the propagating flame front encounters on its way through the containment building. Lumped parameters cannot compute these local turbulence generation processes. Consequently, such codes cannot be used to compute the resulting pressure loads in a reliable way. Furthermore, the propagation of the pressure wave phenomena through the containment cannot be reliably computed using lumped parameter codes due to their inherent limitations (TECDOC, 2011). Therefore, a hydrogen risk assessment methodology based on a combined lumped parameter code/CFD code modeling approach is deemed to be the only viable option to assess the residual risk of possible hydrogen combustion in LWRs. A literature review concerning CFD modeling for hydrogen deflagrations is presented in Sathiah et al. (2012a). NRG’s CFDbased containment model for hydrogen distribution and hydrogen mitigation analyses is described in Houkema et al. (2003) and Visser et al. (2009, 2012a,b, 2013). As explained in Sathiah et al. (2012b), NRG’s CFD-based containment model for hydrogen deflagrations analyses is based on the application of a density based coupled solver for accurate tracking of the induced pressure wave phenomena, the application of an advanced Turbulent Flame speed Closure (TFC) combustion model based on the model of Zimont (1979), and the application of adaptive mesh refinement (AMR) for accurate and efficient tracking of the turbulent flame propagation. In Sathiah et al. (2012b), this modeling approach has been validated against three homogeneous hydrogen–air deflagration experiments performed in the ENACCEF facility (Chaumeix and Bentaib, 2011, 2010; Bleyer et al., 2012; ISP-49, 2011). From the presented validation analyses, it was concluded that the maximum pressures were predicted within 13% accuracy, while the rate of pressure rise dp/dt was predicted within about 30%. The eigen frequencies of the residual pressure wave phenomena were predicted within a few %. Therefore, it was overall concluded that the considered model predicts the considered ENACCEF experiments very well. During a severe accident in a PWR, a hydrogen–air–steam gas mixture will be present in the containment due to the release of both hydrogen and steam from the reactor cooling system. It is important to determine the effect of steam on the hydrogen combustion process. Therefore, we have decided to further extend the validation base of our hydrogen deflagration modeling approach to experiments that include the effect of steam on the hydrogen deflagration process. For that purpose, we have used experimental data from tests that have been executed in the ENACCEF facility. These data have been released within the European SARNET-2 project. In the considered tests, a CO2 –He mixture was used as a substitute diluent in order to avoid condensation effect in the ENACCEF facility. Given this situation, we have defined the objectives of the current paper as follows: • to evaluate the physical properties of hydrogen–air–CO2 –He mixtures and hydrogen–air–H2 O mixtures, in order to assess the effect of using a CO2 –He mixture as a substitute diluent for steam; • to further validate the considered CFD model against three deflagration experiments performed in the ENACCEF facility using different hydrogen–air–CO2 –He mixtures. The paper is structured as follows: Section 2 describes the experiments which have been used for the validation, whereas Section 4 presents the evaluation of the physical properties of hydrogen–air–CO2 –He mixtures and hydrogen–air–H2 O mixtures. Section 3 describes the applied CFD model. The validation results

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

3

2.2. Selected experiments Institut de Radioprotection et de Sureté Nucléaire (IRSN) has released the data of nine ENACCEF experiments for code benchmarking. That is, the data of the tests RUN 153, 158 and 160 (see Table 1) have been released within the SARNET EU project, and the data of three tests with CO2 –He dilution (see Table 1) have been used for code benchmarking within the SARNET-2 project. A 60 vol.% CO2 –40 vol.% He gas mixture was added in these tests to mimic the effect of steam dilution. The data of three tests with non-uniform mixtures (RUN 733, 736 and 765) have been released within the OECD International Standard Problem (ISP-49). In these experiments, the initial temperature is 296 K and the initial pressure is 100 kPa. NRG has used the ENACCEF tests RUN 153, 158 and 160 for the validation of its CFD modeling approach for hydrogen deflagration in Sathiah et al. (2012b). As already mentioned in the introduction, this modeling approach predicted these ENACCEF experiments very well. For the objectives of this paper, three tests with 13 vol.% of hydrogen and respectively 10, 20 and 30 vol.% of CO2 –He added as a substitute diluent for steam (H2 O) have been considered. 3. CFD model

Fig. 1. Schematic representation of ENACCEF geometry.

are presented in Section 5. Finally, the summary and conclusions are presented in Section 6. 2. Selected facility and experiments for code validation 2.1. ENACCEF facility A large number of hydrogen deflagration experiments have been executed in the ENACCEF test facility (Chaumeix and Bentaib, 2011, 2010; Bleyer et al., 2012; ISP-49, 2011). The measurements have been performed by the Centre National de la Recherche Scientifique (CNRS). The test facility consists of a dome and an acceleration tube (see Fig. 1). The acceleration tube is 3.2 m long and has an internal diameter of 0.154 m. It contains 9 annular baffles with varying blockage ratios BR = [1 − (d/D)2 ], where d and D are the baffle and pipe diameter, respectively. The first baffle is located at a distance of 0.776 m and the distance between subsequent baffles is 0.154 m. The thickness of the baffles is 2 mm. The dome is 1.7 m long and has an internal diameter of 0.738 m. In the ENACCEF experiments, the flame position has been measured using 16 UV-sensitive photomultiplier tubes (PMT) (Hamamatsu, 1P28). The presence of the flame was detected based on the total emission of the flame recorded using PMT. Using this method, OH radicals are measured because of their very high concentration in the flame front and in the burnt gases. The uncertainty in the flame position measurements is 8 mm. The pressure signal has been measured using high speed pressure transducers (7 Chimie Metal and 2 Kistler), which are mounted on the inner surface of the tube. Ignition of the gas mixtures in ENACCEF is achieved using two thin tungsten electrodes (2 mm in diameter), which are linked to a high voltage source. The ignition source is located 0.138 m from the bottom of the facility. Further details about the experiments are available in the report of Chaumeix and Bentaib (2011, 2010), Bleyer et al. (2012), and ISP-49 (2011). It is worth stressing here that in the considered experiments, initial turbulence levels were not measured and heat losses from the domain boundary are not specified.

The commercial CFD code ANSYS FLUENT is used as the basis for the development of our hydrogen deflagration model. The unsteady Favre (density-weighted) averaged Navier–Stokes equations are solved using the standard k −  turbulence model. Within the CFD code, the equations for the conservation of mass, momentum, energy and the progress variable are solved. This is already described in our previous paper (Sathiah et al., 2012a). The combustion model with submodels for preferential diffusion, flame stretch effects, and compressions effects on the laminar flame speed for uniform hydrogen–air mixture is described in Sathiah et al. (2012b). Next, the laminar flame speed, the mesh, boundary conditions and the numerics applied in the CFD model are described. 3.1. Laminar flame speed The laminar flame speed of a fuel–air mixture is an important physicochemical property which is needed for combustion modeling. Laminar flame speed measurements for hydrogen–air–diluent mixtures were performed by Liu and MacFarlane (1983) and Bentaib and Chaumeix (2012). Liu and MacFarlane (1983) measured the laminar flame speed of hydrogen–air–steam mixtures in a Bunsen burner. The experiment covers hydrogen concentrations in the range of 18–65% (by vol.), H2 O concentrations in the range of 0–15% (by vol.), a temperature range of 296–523 K, and a pressure of 1 atm. They fitted their measured laminar flame speed into an empirical formula as follows: Sl,0 = BT C exp(DX steam ),

(1)

where, B, C and D are constants given by B = A1 + A2 (0.42 − XH2 ) + A3 (0.42 − XH2 )2 ,

(2)

and C = A4 + A5 (0.42 − XH2 ) and D = A6 ,

(3)

where Xsteam is the mole fraction of H2 O, XH2 is the mole fraction of hydrogen, T (K) is the temperature of the unburned gas, and A1 , A2 , A3 , A4 , A5 and A6 are constants. The values of the constant A1 . . . A6 are summarized in Liu and MacFarlane (1983). Bentaib and Chaumeix (2012) measured the laminar flame speed for hydrogen–air–CO2 –He mixtures for different diluent

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model

ARTICLE IN PRESS

NED-7927; No. of Pages 15

P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

4

Table 1 ENACCEF experiments released for code benchmarking. RUN

[H2 ] vol.%

Blockage ratio

Gas mixture

153 158 160 153-10 153-20 153-30 765 736 733

13 13 13 13 13 13 11.6–8.1 11.4–5.8 5.7–12

0.63 0.33 0 0.63 0.63 0.63 0.63 0.63 0.63

Uniform hydrogen–air mixture Uniform hydrogen–air mixture Uniform hydrogen–air mixture Uniform hydrogen–air mixture with 10 vol.% CO2 –He Uniform hydrogen–air mixture with 20 vol.% CO2 –He Uniform hydrogen–air mixture with 30 vol.% CO2 –He Hydrogen–air mixture with negative concentration gradient Hydrogen–air mixture with negative concentration gradient Hydrogen–air mixture with positive concentration gradient

mole fractions. They claimed that a mixture of 60% CO2 and 40% He mixture can mimic H2 O as diluent, since both of these diluents have similar physical properties (Lamoureux et al., 2003). The measurement results are fitted into the following correlation for the laminar flame speed:



Sl,0 = (1.442 + 1.07 − 0.29)(1 −

)

4

T Tref

2.2 

P Pref

−0.5

, (4)

is the mole fraction of the where  is the equivalence ratio, diluent, Pref is the reference pressure and Tref is the reference temperature which take values of 100 kPa and 298 K respectively. In the above correlation, the pressure P varies between 100 and 500 kPa, while the temperature T varies between 298 and 353 K. Fig. 2 shows the variation of the laminar flame speed with diluent concentration as calculated using the two correlations given in Eqs. (1)–(4). It can be observed that an increase in the diluent mole fraction decreases the laminar flame speed. The laminar flame speed values obtained by Liu and MacFarlane (1983) are higher than obtained using the correlation proposed by Bentaib and Chaumeix (2012). This is because the correlation proposed by Liu and MacFarlane (1983) is not corrected for stretch effects (Poinsot and Veynante, 2001). For our simulations presented in Section 5, we extracted the laminar flame speed values from the correlation proposed by Bentaib and Chaumeix (2012), i.e. Eq. (4). Moreover, the effects of compression of the unburned gases on laminar flame speed must be taken into account as this important for explosion modeling. This is included in Eq. (4). 3.2. Mesh, initial and boundary conditions The ENACCEF facility is represented using a 2D axisymmetric geometry. As described in Sathiah et al. (2012b), an AMR method is used to resolve the flame brush thickness. The base grid used in the simulations is the same as the one used in our earlier work (Sathiah et al., 2012b). The walls of the domain were assumed to

be adiabatic, and no flame-wall interaction was taken into account. The mesh dimensions in x and y directions before and after AMR are 0.02 m and 0.0025 m, respectively. Detailed grid and time step sensitivity analyses have been performed for each case in order to guarantee that the numerical errors can be practically neglected in the present validation analyses. Stagnant flow conditions have been applied as initial conditions. Since the initial turbulence parameters have not been measured in the ENACCEF experiments, best guess values of 1.5e−04 m2 s−2 and 4.8e−05 m2 s−3 have been assumed for the turbulent kinetic energy and the turbulent dissipate rate respectively. The corresponding values of the turbulent length scale and turbulent intensities are 0.007 m and 0.01 m/s, respectively. These values were also used for the three experiments considered in our previous paper (Sathiah et al., 2012b). It is worth stressing that our best guess values for initial turbulence are not tuned to match the experimental results.

3.3. Models and postprocessing method Since the ENACCEF facility is a vertical facility, buoyancy effects are important. Effects of buoyancy on turbulent kinetic energy and turbulent dissipation rate were accordingly considered in the simulation. We used an ignition radius of 1 cm in the present analyses. As demonstrated in Sathiah et al. (2012a), a small increase in ignition radius results only in a small time shift of the results. In addition, it is assumed in our simulations that counter-gradient transport has negligible effects on the results. Effects of stretching on flame propagation were not considered, which means that the function G in Eq. (4) in Sathiah et al. (2012b) is set to unity. However, preferential diffusion and compression effects were considered in the simulations. The flame position obtained from the simulations was extracted from the axial-coordinate corresponding to an iso-surface value c˜ = 0.5, as done by Gubba et al. (2009). 3.4. Applied numerical schemes

Fig. 2. The variation of the laminar flame speed Sl,0 with diluent concentrations measured by Bentaib and Chaumeix (2012) and Liu and MacFarlane (1983) for a hydrogen concentration of 13 vol.%.

The density based coupled solver of ANSYS FLUENT (2008) is used in our simulations. It solves the conservation equations for mass, momentum, energy, and the progress variable (Eqs. (1)–(4) in Sathiah et al., 2012a) simultaneously. The equations for the turbulent kinetic energy and the turbulent dissipation rate (Eqs. (7) and (8) in Sathiah et al., 2012a) are solved sequentially. The code employs a 4-stage Runge–Kutta scheme which is fourth order accurate for time integration of unsteady flows. The time step in this method is determined by the Courant–Friedrichs–Lewy (CFL) condition. This time-step is used in each control volume of the domain. A CFL number of 0.8 is used in the simulations. A complete description is provided in Blazek (2005). The spatial discretization of flow, species and turbulence equations is performed using the second-order upwind numerical scheme. The convective flux is calculated using the Flux-Difference Splitting (FDS) scheme by Roe (1981, 1986).

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

Fig. 3. The variation of the unburnt molecular diffusivity with diluent concentration (vol.%) calculated for hydrogen–air mixtures with the different diluents H2 O and CO2 –He.

5

Fig. 5. The variation of the unburnt gas Lewis number Le with diluent concentration (vol.%) calculated for hydrogen–air mixtures with the different diluents H2 O and CO2 –He.

This trend is caused by the fact that the thermal conductivity of CO2 –He is larger than the thermal conductivity of H2 O. This different behavior of diluents can also be observed for the Lewis number (see Fig. 5). A 40% increase in the diluent mole fraction increases the Lewis number by 16% for a hydrogen–air–CO2 –He mixture, while it decreases by 8% for a hydrogen–air–H2 O mixture. Figs. 6 and 7 demonstrate the variation of the molecular weight of the unburnt and burnt gases with the diluent mole fraction. The effects of an increase in the diluent mole fraction on unburnt and burnt molecular weight is significant for hydrogen–air–H2 O mixtures. This is because of the lower molecular weight of H2 O (18 g/mol) in comparison to CO2 –He (44 g/mol). Because of this lower molecular weight, the unburnt gas density Fig. 4. The variation of the unburnt thermal diffusivity u with diluent concentration (vol.%) calculated for hydrogen–air mixtures with the different diluents H2 O and CO2 –He.

4. Comparison of the physical properties of hydrogen–air–H2 O and hydrogen–air–CO2 –He mixtures As mentioned in the previous section, a mixture of 60% CO2 and 40% He (CO2 –He) was used as a diluent instead of steam (H2 O) in the ENACCEF experiments. In order to assess the effects of using a CO2 –He mixture instead of H2 O as diluent, it is necessary to compare the physical properties of hydrogen–air mixtures containing these diluents. The physical properties compared here are divided into the ones which affects (a) the flame dynamics, (b) the pressure dynamics, and (c) the acoustic phenomena observed in the ENACCEF experiments. The physical properties compared here are calculated using GASEQ (2011) and CANTERA (2011) tools. The comparative analyses are presented below.

Fig. 6. The variation of the unburnt molecular weight with diluent concentration (vol.%) calculated for hydrogen–air mixtures with the different diluents H2 O and CO2 –He.

4.1. Physical properties affecting the flame dynamics Fig. 3 illustrates the variation of the unburnt molecular diffusivity (i.e. diffusivity of hydrogen in the unburnt mixture) with the diluent mole fraction. It can be noticed that the molecular diffusivities of hydrogen–air mixtures with the two considered diluents are quite similar. The molecular diffusivity increases with increase in the diluent mole fraction. A 40% increase in diluent mole fraction increases the molecular diffusivity by 6%. The corresponding thermal diffusivity shows a different behavior. Namely, it increases with increase in the mole fraction for hydrogen–air–CO2 –He mixtures, while for hydrogen–air–H2 O mixtures, it decreases with increasing diluent concentration (see Fig. 4). This opposite trend is due to the increase in the thermal conductivity of the mixture with increase in the diluent mole fraction for hydrogen–air–CO2 –He mixtures.

Fig. 7. The variation of the burnt molecular weight with diluent concentration (vol.%) calculated for hydrogen–air mixtures with the different diluents H2 O and CO2 –He.

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15 6

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

Fig. 8. The variation of the unburnt gas density u with diluent concentration (vol.%) calculated for hydrogen–air mixtures with the different diluents H2 O and CO2 –He.

Fig. 10. The variation of the absorption coefficients of the burnt gases as a function diluent mole fraction as calculated for hydrogen–air mixtures with two different diluents H2 O and CO2 –He.

strongly decreases with increase in the diluent mole fraction in case of hydrogen–air–H2 O mixtures (see Fig. 8). That is a 40% increase in the diluent mole fraction decreases the unburnt gas density by 20% for a hydrogen–air–H2 O mixture, while it decreases by 0.1% for a hydrogen–air–CO2 –He mixture. Fig. 9 shows the variation of the specific heat capacity of hydrogen–air–CO2 –He and hydrogen–air–H2 O mixtures. The increase in the specific heat capacity of the mixture with the diluent mole fraction is stronger for hydrogen–air–H2 O mixtures than for hydrogen–air–CO2 –He mixtures. This is because of the higher specific heat capacity of H2 O in comparison to CO2 –He. That is, the specific heat capacity increases by about 30% for 40% H2 O dilution and by about 3% for 40% CO2 –He dilution. The heat loss by thermal radiation is obtained from the absorption coefficient of the burnt gas. Fig. 10 compares the absorption coefficient of the burnt gases containing two different diluents as a function of diluent mole fraction. The absorption coefficient is calculated at adiabatic flame temperature. The data for the absorption coefficient is taken from Beyler et al. (2002). It is worth mentioning that the absorption coefficient of He is zero. It can be observed that the absorption coefficient of the burnt gases containing CO2 –He diluent is smaller than the absorption coefficient of the burnt gases containing H2 O as diluent. Additionally, the difference in the absorption coefficient increases with increase in dilution. A 30% increase in dilution increases the absorption coefficient of the burnt gases by 276% for CO2 –He diluent and by 322% for H2 O diluent. This essentially means that the radiation effects of hydrogen–air–H2 O and hydrogen–air–CO2 –He mixtures are different. To conclude, differences exist in the molecular, thermal diffusivity and Lewis number, unburnt gas density, specific heat capacity and absorption coefficient of the gases for hydrogen–air–CO2 –He

The above mentioned properties are investigated because it can affect the flame dynamics, i.e. the laminar and turbulent flame propagation. Furthermore, the effects of using different diluents on the pressure dynamics, i.e. on the rate of pressure rise and on the maximum pressure must be investigated. For this purpose, the adiabatic flame temperature and AICC (Adiabatic Isochoric Complete Combustion) pressure obtained for hydrogen–air–H2 O and hydrogen–air–CO2 –He mixtures as a function of diluent mole fractions are compared. Figs. 11 and 12 show the variation of the adiabatic flame temperature and the AICC pressure with diluent concentration for hydrogen–air–H2 O and hydrogen–air–CO2 –He mixtures. It can be noticed that the adiabatic flame temperature decreases with increase in the diluent concentration. This is caused by the increase in the specific heat capacity of the mixture with an increase in the diluent mole fraction (see Fig. 9). A 40% increase in diluent mole fraction decreases the adiabatic flame temperature by 6.8% for CO2 –He, while it decreases by 8.3% for H2 O. This also influences the AICC pressure, which also decreases with increase in the diluent mole fraction. The considered decrease in the adiabatic flame temperature and AICC pressure are slightly larger for the H2 O diluent. Namely, a 40% increase in diluent mole fraction decreases the AICC pressure by 6.8% for CO2 –He, while it decreases by 8.3% for H2 O. This means that the values of maximum temperature

Fig. 9. The variation of the specific heat capacity Cp as a function of diluent mole fraction calculated for hydrogen–air mixtures with the different diluents H2 O and CO2 –He.

Fig. 11. The variation of the adiabatic flame temperature Tad with diluent concentration (vol.%) calculated for hydrogen–air mixtures with the different diluents H2 O and CO2 –He.

and hydrogen–air–H2 O. This essentially means that the laminar and turbulent flame propagation are different in both these mixtures. 4.2. Physical properties affecting the pressure dynamics

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

Fig. 12. The variation of the AICC pressure pAICC with diluent concentration (vol.%) calculated for hydrogen–air mixtures with the different diluents H2 O and CO2 –He.

and pressure are slightly larger in hydrogen–air–CO2 –He mixtures than in hydrogen–air–H2 O mixtures.

7

Fig. 13. The variation of the speed of sound c with diluent concentration (vol.%) calculated for hydrogen–air mixtures with the different diluents H2 O and CO2 –He.

hydrogen–air–H2 O mixtures, especially at high diluent concentrations, provided the acoustic intensities are the same. 4.4. Evaluation

4.3. Physical properties affecting the pressure wave acoustic phenomena

The dilution of hydrogen–air experiments with CO2 –He results in:

In addition, for the pressure wave acoustic phenomena observed in the ENACCEF experiments, the speed of sound, and the acoustic impedance are of importance, since these two parameters determine the frequencies and amplitudes of the residual pressure waves after complete burning. For example, let us consider the case of acoustic propagation in a pipe. The eigen frequency ω of the pipe is a function of the speed of sound, namely: ω=

c , 4L

(5)

where c is the speed of sound, and L is the length of the pipe. Recently, this formula was also used by Movahed-Shariat-Panahi (2012) to estimate the mesh size needed to resolve pressure oscillations. Furthermore, the amplitude of the observed pressure waves (ıp) is proportional to the acoustic impedance (Z) (Poinsot and Veynante, 2001) as follows ıp ∝ Z = b c, ıu

(6)

where ıu is the velocity fluctuation and b is the burnt gas density. It can be observed from the above expression that an increase in acoustic impedance increases the amplitude of pressure fluctuations provided that the acoustic intensity I = ıpıu is constant. Therefore, it is worth comparing the speed of sound and the acoustic impedance of hydrogen–air–CO2 –He and hydrogen–air–H2 O mixtures. Fig. 13 shows an increase in the speed of sound with an increase in the diluent mole fraction for hydrogen–air–H2 O mixtures, while the speed of sound decreases for hydrogen–air–CO2 –He mixtures. The increase in the speed of sound in case of hydrogen–air–H2 O mixtures is due to the strong decrease in the molecular weight of the mixture with an increase in the diluent mole fraction (see Fig. 7). A 40% increase in diluent mole fraction decreases the speed of sound by 4% for CO2 –He diluent, while it increases the speed of sound by 3.5% for H2 O diluent. Fig. 14 shows a decrease in the acoustic impedance with an increase in the diluent concentration. The decrease in the acoustic impedance for hydrogen–air–CO2 –He mixtures is smaller than for hydrogen–air–H2 O mixtures. A 40% increase in diluent mole fraction decreases the acoustic impedances by 5% and 14.5% for the CO2 –He and H2 O diluents respectively. This essentially means that the amplitude of the pressure fluctuations for hydrogen–air–CO2 –He mixtures are larger than for

• • • •

a lower laminar flame speed, a lower adiabatic flame temperature, a lower AICC pressure, smaller eigen frequencies and amplitudes of the residual pressure oscillations. More specifically, a 30% dilution of CO2 –He results in

• a decrease in the laminar flame speed by about 37%, • a decrease in the adiabatic flame temperature and AICC pressure by about 5.1%, • a decrease in the eigen frequencies and amplitude of the pressure oscillations by 3%. while in case of H2 O diluent a 30% dilution results in • a decrease in the laminar flame speed by about 48.5%, • a decrease in adiabatic flame temperature and AICC pressure by 6.2%, • and increase in the eigen frequencies by 2.6% and a decrease in the amplitude of the pressure oscillations by 10.8%. To summarize, the physical properties, the adiabatic flame temperature, the AICC pressure, the speed of sound, and the

Fig. 14. The variation of the acoustic impedance Z = b c with diluent concentration (vol.%) calculated for hydrogen–air mixtures with the different diluents H2 O and CO2 –He.

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

8

acoustic impedance of hydrogen–air mixtures with CO2 –He and H2 O as diluents are somewhat different. Therefore, the flame and the pressure dynamics in hydrogen–air mixtures with CO2 –He and H2 O as diluents will be somewhat different. This essentially means that experimentalist should be careful in using CO2 –He as a diluent to mimic the effect of steam diluent on the flame and pressure dynamics of hydrogen–air mixtures. Furthermore, the modelers should be careful in using the experimental data of hydrogen–air–CO2 –He for the validation of the code to simulate combustion of a hydrogen–air–H2 O mixtures. To check the capability of our model to show the differences in the flame and pressure dynamics between hydrogen–air–CO2 –He and hydrogen–air–H2 O mixtures, simulations have been performed for three different diluent mole fractions. These results are presented in Section 5. 5. Validation results The CFD model presented in Section 3 is validated against deflagration experiments performed in the ENACCEF facility for a uniform hydrogen–air mixture with different diluent concentrations. For this purpose, we selected four different cases namely RUN 153, 153-10, 153-20 and 153-30 (see Table 1). The data of these tests were released in the framework of the European SARNET-2 project. In these experiments, a mixture of CO2 and He is applied as a diluent to substitute H2 O in order to avoid possible condensation effects. In Section 4, it was shown that dilution by steam has somewhat bigger effects on the AICC pressure and adiabatic flame temperature compared to the dilution with CO2 –He. Therefore, computations have been performed also for steam dilution, in order to determine whether our modeling predicts these trends qualitatively well. In the following subsections, the simulation and experimental results are compared for RUN 153-10, RUN 153-20 and RUN 15330. Because of the importance of grid and time-step on the flame propagation, detailed sensitivity studies were performed for each case to guarantee grid and time-step independent results. For completeness, the validation results corresponding to 0% dilution (i.e. RUN 153 in Table 1) are also presented here. The corresponding results were already presented in our previous work (Sathiah et al., 2012b). 5.1. RUN 153 Fig. 15a shows the variation of the axial flame position with time obtained for RUN 153 without diluent. The numerical results are shown for three different levels of AMR. It can be concluded that practically identical results were obtained using 2 and 3 levels of AMR. This can be further confirmed from the computed pressures (see Fig. 15b). Hence, it can be concluded that 2 levels of AMR is sufficient to obtain practically grid independent results. The grid independence was further confirmed by comparing the volume integral source term (right hand side of progress variable equation) for three different levels of AMR. However, these results are not shown here for brevity reasons. The following four different phases of flame propagation (see Fig. 16a) can be identified in ENACCEF: (a) the initial quasi-laminar flame propagation, (b) turbulent flame acceleration, (c) flame deceleration, and (d) the jet flame. The first phase starts when the flame is ignited and lasts until the flame reaches the first baffle, more specifically at time t = 0–0.025 s in the simulations. The times at which the flame reaches respectively the first baffle, the last baffle, and the dome are also shown here. In the quasi-laminar flame propagation phase, the flame propagates slowly because of low turbulence levels. The flame velocity (defined as the ratio of axial flame

Fig. 15. Grid sensitivity results for RUN 153 (0 vol.% dilution% and blockage ratio = 0.63). The variation of (a) the axial flame front position and (b) the pressure with time calculated using three different levels of AMR.

position s with time t, i.e. ds/dt) predicted by the simulation in this phase is higher than in the experiments. The possible causes for the overprediction could be (a) even lower initial turbulence levels in the experiments than that we have assumed. (b) overprediction of the turbulence levels by the standard k −  turbulent model in the quasi-laminar regime, (c) overprediction of the turbulent source term by the applied TFC combustion model in the quasi-laminar regime. This overprediction in the flame velocity during the initial quasi-laminar flame propagation phase essentially introduces a time shift between simulation results and experiments. As will be seen later, this overprediction does not affect the turbulent flame acceleration phase and corresponding pressure dynamics. The subsequent turbulent flame acceleration phase starts when the flame reaches the first baffle and ends when the flame arrives at the last baffle in the acceleration tube (i.e. time t = 0.025–0.036 s). The flame velocity in this phase increases rapidly because of the turbulence generated by the baffles. The mechanisms causing the increased burning rate are the wrinkling of the flame front by the generated turbulent eddies and the turbulent transport of heat and mass at the reaction front. The simulation very well predicts this phase and the slope of the axial flame front position with time compares well with the experiments during this turbulent flame acceleration phase. The third phase is the flame deceleration phase which starts when the flame passes the final obstacle and ends when it reaches the dome entrance (i.e. time t = 0.036–0.042 s). In this phase, the flame velocity decreases due to a decrease in the turbulence levels. The slope of the flame distance versus time as predicted by simulations is in good agreement with experiments.

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

9

Table 2 Quantitative comparison of experimental and numerical data for ENACCEF experiments with 0, 10, 20, and 30% of dilution. + and − indicate respectively the overprediction and underprediction with respect to the experimental value.

RUN 153 First peak pressure Maximum pressure First rate of pressure rise dp/dt Second rate of pressure rise dp/dt First eigen frequency Second eigen frequency Amplitude corresponding to first eigen freq. Amplitude corresponding to second eigen freq. RUN 153-10 First peak pressure Maximum pressure First rate of pressure rise dp/dt Second rate of pressure rise dp/dt First eigen frequency Second eigen frequency Amplitude corresponding to first eigen freq. Amplitude corresponding to second eigen freq. RUN 153-20 First peak pressure Maximum pressure First rate of pressure rise dp/dt Second rate of pressure rise dp/dt First eigen frequency Second eigen frequency Amplitude corresponding to first eigen freq. Amplitude corresponding to second eigen freq.

Fig. 16. Results for RUN 153 (0 vol.% dilution and blockage ratio = 0.63). The variation of (a) the axial flame front position and (b) the pressure with time and the variation of (c) the amplitude with frequency calculated using 2 levels of AMR.

The final phase is the turbulent flame expansion where the flame resembles a jet. The flame radially expands into the dome, pushes the unburnt gases present in the dome and slowly consumes the remaining unburnt gases. During the quasi-laminar and turbulent flame acceleration phases, the rise in the pressure is negligible (see Fig. 16b). During the flame acceleration, compression pressure waves detached from the flame and travels to the dome entrance. At this dome entrance the compression waves are reflected as expansion waves. The reflected expansion waves create the sudden decrease in the intermediate peak pressure during the deceleration phase. The intermediate peak in the pressure is observed at 0.0375 s in the simulation, while in the experiments this is observed at 0.9 s. As indicated, this time shift is caused by the faster flame propagation during the quasi-laminar flame propagation phase in the

RUN 153-30 First peak pressure Maximum pressure First rate of pressure rise dp/dt Second rate of pressure rise dp/dt First eigen frequency Second eigen frequency Amplitude corresponding to first eigen freq. Amplitude corresponding to second eigen freq.

Experiment

TFC

Difference

2.49e05 Pa 4.5643e05 Pa 3.5929e09 Pa/s

2.62e05 Pa 4.7638e05 Pa 3.01679e09 Pa/s

+5.2% +4.3% −16%

1.72e07 Pa/s

1.9961e07 Pa/s

226.7 Hz 460.3 Hz

231.3 Hz 456.3 Hz

5590 Pa

2679 Pa

−51.9%

2114 Pa

1329 Pa

−37.3%

3.562e05 Pa 4.353e05 Pa 1.86e10 Pa/s

2.5255e05 Pa 4.868e05 Pa 2.071e09 Pa/s

−29.1% +11.8% −88.8%

2.0543e07 Pa/s

1.7543e07 Pa/s

−14.6%

228.6 Hz 457.1 Hz

232.8 Hz 456.9 Hz

4856 Pa

5367 Pa

4443 Pa

887 Pa

3.357e05 Pa 4.339e05 Pa 1.537e10 Pa/s

2.394e05 Pa 4.67e05 Pa 1.721e09 Pa/s

−28.6% +7.6% −88.8%

2.321e07 Pa/s

1.5041e07Pa/s

−35.1%

228.6 Hz 457.1 Hz

230.6 Hz 471.3 Hz

4956 Pa

2201 Pa

−55.5%

8378 Pa

189 Pa

−97.7%

2.73e05 Pa 4.253e05 Pa 5.99e09 Pa/s

2.115e05 Pa 4.492e05 Pa 8.2815e08 Pa/s

−22.5% +5.6% −86.1%

2.76e07 Pa/s

1.173e07 Pa/s

−57.5%

220 Hz 446.7 Hz

228.1 Hz 473.3 Hz

9276 Pa

1808 Pa

−80.5%

3750 Pa

293.9 Pa

−92.1%

+15.7% +2.0% −0.8%

+1.8% −0.04% +10.5%

−80%

+0.8% +3.1%

+3.6% +5.9%

simulation which introduces a time-shift of approximately 0.6 s. Furthermore, the values of the first intermediate peak pressure is slightly overpredicted in the simulation in comparison to the experiments (see Table 2).

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15 10

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

Followed by the first intermediate peak in the pressure, the pressure increases rapidly because of the burning of the gases in the dome, until it reaches the maximum pressure. This increase in the pressure with time, i.e. the slope of the pressure dp/dt is captured well. The maximum value of the mean pressure is slightly overpredicted, because adiabatic boundary condition is used in the simulations. Once the maximum pressure is reached, the pressure decreases with time in the experiments because of heat loss from the facility. Because heat loss is not incorporated in the applied CFD model, this slow pressure decrease is not captured in the simulations. Once the pressure has reached its maximum value, significant residual pressure oscillations occur due to the residual pressure wave phenomena observed in the facility. Such pressure oscillations can potentially jeopardize the structural integrity of the containment or safety systems inside the containment when the frequencies of the pressure oscillations match the eigen frequencies of the containment or safety systems. To analyze the frequency spectrum of the pressure, Fast Fourier Transform (FFT) analyses were performed. Fig. 16c shows the variation of the amplitude with frequency obtained using simulations and experiments. It is worth mentioning that since the amplitude of the oscillations decreases with time, the values of the amplitude shown here are time averaged. The first and second frequency are predicted very well by the simulations while the corresponding amplitudes are slightly underpredicted (see Table 2). 5.2. RUN 153-10 Fig. 17a shows the variation of the axial flame front position with time obtained using three different levels of AMR. For brevity, only the results for the hydrogen–air–CO2 –He mixture are presented here. Again, practically grid independent solutions were obtained by using 2 levels of AMR. The variation of the pressure with time also shows similar grid independence for 2 levels of AMR (see Fig. 17b). The same conclusion was drawn for the hydrogen–air–H2 O mixture. It is worth mentioning that the base grid used in the present simulation is the same as the one used in our previous ENACCEF simulation (see Sathiah et al., 2012b). There we obtained also grid independent solutions with two levels of AMR. Furthermore, it was concluded that further 2 by 2 refinement of the original base grid had no effect on flame and pressure dynamics. Hence, this study was not repeated in the present work. The variation of the axial flame front position with time is shown in Fig. 18a for the different diluents. In the experiments, it is assumed that the hydrogen–air mixture was ignited at time t = 0 s. As for RUN 153 (see Section 5.1), the same four phases of flame propagation can be identified. The flame velocity in the simulations is again overpredicted during the quasi-laminar flame propagation phase. This overprediction introduces a time shift of 0.06 s between the experiment and simulations. Similar as for RUN 153, this timeshift is due to overprediction in the flame velocity during the initial quasi-laminar phase. However, the turbulent flame acceleration and flame deceleration phases are well predicted. The predicted values of the intermediate pressure peak is lower than in the experiment (Fig. 18b). These values are quantified in Table 2. During the flame deceleration phase, the pressure increases quickly with time because of the combustion of the unburned gases present in the dome. This phase is predicted well in comparison to the experiments. The corresponding values of the first and second dp/dt as predicted by the simulations are provided in Table 2. The values of the maximum mean pressure obtained by the simulations is somewhat higher than in the experiments. Again, this is because of the heat loss effects which were not considered in the simulations, since adiabatic boundary conditions were used. As a result, the decrease in the pressure with time as observed in

Fig. 17. Grid sensitivity results for RUN 153-10 (10 vol.% dilution and blockage ratio = 0.63). The variation of (a) the axial flame front position and (b) the pressure with time calculated using three different levels of AMR. The results correspond to the hydrogen–air–CO2 –He mixture.

the experiments is not reproduced in the simulations. However according to the authors’ opinion, the prediction of this decrease in pressure with time is not important for the safety analyses of the containments. A lower value of maximum value of the mean pressure is also observed in the simulations corresponding to the hydrogen–air–H2 O mixture. The reason for this is due the larger value of specific heat capacity of H2 O diluent, as described in Section 4 (see Fig. 12). Fig. 18c shows the variation of the amplitude with frequency. The first and second frequency predicted by the simulations compares very well with the experiments. However, the model slightly overpredicts the amplitudes corresponding to the first eigen frequency and strongly underpredicts the amplitude of the second eigen frequency (see Table 2). The frequency and corresponding amplitudes predicted by the simulations for hydrogen–air–H2 O and hydrogen–air–CO2 –He mixture are very similar. 5.3. RUN 153-20 The results of the grid sensitivity performed for RUN 153-20 are shown in Fig. 19a and b. From the presented results, it was concluded that practical grid independent results were obtained for two levels of AMR. Fig. 20a shows the variation of axial flame front with time for the simulations with the different diluents. During the quasi-laminar phase, the flame velocity is again overpredicted in the simulations. The reason for this is already explained before. This eventually leads to a time shift of 0.067 s between the experiments and simulations. However, the simulations predict the rapid increase in the

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

11

Fig. 19. Grid sensitivity results for RUN 153-20 (20 vol.% dilution and blockage ratio = 0.63). The variation of (a) the axial flame front position and (b) the pressure with time calculated using three different levels of AMR for the hydrogen–air–CO2 –He mixture.

first eigen frequency. This trend is not observed in the RUN 153-10 and RUN 153-30. The reason for this behavior in RUN 153-20 is not clear to the authors. The differences in the frequencies and corresponding amplitudes between the different diluents CO2 –He and H2 O are small. 5.4. RUN 153-30 Fig. 18. Results for RUN 153-10 (10 vol.% dilution and blockage ratio = 0.63). The variation of (a) the axial flame front position, (b) the pressure with time and the variation of (c) the amplitude with frequency calculated using 2 levels of AMR for hydrogen–air–CO2 –He and hydrogen–air–H2 O mixtures.

turbulent flame velocity due to the turbulence generated during the turbulent flame acceleration phase well. The intermediate peak pressures obtained in the simulations are somewhat underpredicted in comparison to the experiment (Fig. 20b). The predicted intermediate peak pressure observed for the CO2 –He and H2 O diluents are very similar. Furthermore, the maximum value of the mean pressure reached in the simulations is slightly higher than in the experiments. The maximum mean pressure for the hydrogen–air–CO2 –He mixture is slightly higher than the maximum mean pressure for hydrogen–air–H2 O mixture. The reason for this difference is described in Section 4 (see Fig. 12). The pressure spectrum in Fig. 20c shows that the first and second frequency are predicted very well. The amplitude corresponding to the first eigen frequency is underpredicted in the simulations, while it is strongly underpredicted for the second eigen frequency. It can be observed in the experiments that the amplitude corresponding to the second eigen frequency is higher than the amplitude for the

Fig. 21a and b illustrate the variation of the axial flame front position and the pressure with time obtained using three different levels of AMR. Strictly speaking, it cannot be concluded that grid independent results are obtained. However basically, there is just a very small time shift between the results. Therefore similar as in the previous cases, we have based our subsequent analyses on the results obtain with 2 levels of AMR. Fig. 22a demonstrates the variation of the axial flame front position with time obtained for a diluent concentration of 30 vol.%. As in the case of RUN 153, 153-10 and 153-20, the overprediction of the flame velocity during the quasi-laminar phase is also observed in this case. This leads to time-shift of 0.095 s between the experiment and simulations. Fig. 22b shows the variation of the pressure with time. As in the case of RUN 153, 153-10 and 153-20, the intermediate peak pressure is underpredicted in the simulations. Next, the pressure increases quickly with time when the flame is present in the dome. It is observed that the maximum value of the mean pressure is also overpredicted in the simulations due to the application of adiabatic boundary conditions. This is consistent with the results of RUN 153, 153-10 and 153-20.

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15 12

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

Fig. 20. Results for RUN 153-20 (20 vol.% diluent and blockage ratio = 0.63). The variation of (a) the axial flame front position with time, (b) the pressure with time and (c) the amplitude with frequency calculated using 2 levels of AMR for hydrogen–air–CO2 –He and hydrogen–air–H2 O mixtures.

Fig. 22. Results for RUN 153-30 (30 vol.% dilution and blockage ratio = 0.63). The variation of (a) the axial flame front position, (b) the pressure with time and the variation of (c) the amplitude with frequency calculated using 2 levels of AMR for hydrogen–air–CO2 –He and hydrogen–air–H2 O mixtures.

Fig. 21. Grid sensitivity results for RUN 153-30 (30 vol.% dilution and blockage ratio = 0.63). The variation of (a) the axial flame front position and (b) the pressure with time calculated using three different levels of AMR for the hydrogen–air–CO2 –He mixture.

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

13

Fig. 22c shows the variation of the amplitude with frequency. The first and second eigen frequency are predicted well. However, their corresponding amplitudes are strongly underpredicted (see also Table 2). In the above results, it can be noticed that the flame and pressure dynamics as predicted by the simulations for the H2 O and CO2 –He as diluents are very similar. Somewhat larger differences occur for the maximum value of the mean pressure. 5.5. Qualitative comparison In order to analyze the effects of dilution on flame dynamics, pressure dynamics, and pressure wave phenomena in a qualitative way, the experimental results and simulations results of RUN 153, 153-10, 153-20, and 153-30 are presented here group wise. Fig. 23 shows axial flame front position with time measured in the experiments for four different diluent concentrations. The experimental results of RUN 153 show inconsistent behavior in comparison to other experiments. This experiment was a part of a different batch of experiments. As a result, the observed difference could be due to a different ignition time delay in igniting the unburnt gases. It can be observed from the slope of the axial flame position versus time that the flame velocity decreases with an increase in the diluent mole fraction (10–30%), as shown in Fig. 23a. The flame velocity decreases due to a decrease in the laminar flame speed with an increase in the diluent mole fraction (see Fig. 2). This trend is very well captured in the simulations (see Fig. 24a). The slope of the axial flame position versus time predicted by the simulations for RUN 153, 153-10, 153-20 and 153-30 are 272.3, 260.8, 204.3 and 183 m/s, respectively. For experiments, the slope of the axial flame position for RUN 153, 153-10, 153-20 and 153-30 are 386.44, 379.23, 328.08 and 315.9 m/s, respectively. This means that the simulations are able to predict the trend of decrease in the flame velocity with an increase in the dilution concentration well. In the experiments, the intermediate peak pressure decreases with an increase in the diluent mole fraction (10–30%) as shown in Fig. 23b. This trend is also very well reproduced in the simulations (see Fig. 24b). The experiments show a decrease in the first dp/dt when the dilution concentration is increased from 10 to 30% (in contrast, the first dp/dt increases when dilution is increased from 0 to 10%). The CFD simulations also show a decrease in the first dp/dt with increase in the diluent fraction. This decrease in the first dp/dt is due to the decrease in the laminar flame speed with an increase in diluent mole fraction. In contrast, the second dp/dt increases in the experiments, while the CFD simulations show a decrease in the second dp/dt. So, the CFD model does not capture this trend. This opposite trends takes place during the fourth phase of the flame propagation, that is during the turbulent jet flame in the dome. Since the experimental data for the flame position versus time is missing, we were not able to make a further investigation of this fourth phase. In the experiments, the maximum value of the mean pressure decreases with an increase in the diluent mole fraction. This is due to the larger heat capacity of the diluent which decrease the adiabatic flame temperature, thereby, reducing the maximum mean pressure (see Fig. 12). This trend is captured well in the simulations. This means that simulations are able to predict the trend of increase in the dilution on the intermediate peak pressure, the first dp/dt, and the maximum value of the mean pressure as observed in the experiments. It is worth mentioning that these trends are consistent with the trends presented in Section 4. In particular, see Fig. 12, which shows the decrease in the AICC pressure with an increase in dilution. It can be noticed in the experiments that the increase in the diluent mole fraction has very little effect on the values of the first and second eigen frequency as shown in Fig. 23c. This trend is also

Fig. 23. The variation of (a) the axial flame front position, (b) the pressure with time and the variation of (c) the amplitude with frequency obtained using experiments for hydrogen–air–CO2 –He for four different diluent concentrations. The results corresponds to experiments.

observed in the simulations, as shown in Fig. 24c. Moreover, these trends are consistent with trends presented in Section 4, where we observed that dilution effects the speed of sound by just few %. In the experiments, the amplitudes corresponding to the first eigen frequency, first decrease and then increase with an increase in the diluent mole fraction. In contrast, the model does not reproduce this trend. That is the model predicts a systematic decrease in the amplitude of these oscillations with an increase in diluent mole fraction, which is consistent with the trends presented in Section 4. Namely, according to the theory presented in Section 4, an increase in the diluent fraction results in a decrease in the acoustic impedance Z (see Fig. 14) which in turn results in smaller amplitude of the pressure oscillations.

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15 14

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

parameters for the four different cases RUN 153, 153-10, 153-20, and 153-30. From this table, it can be concluded that the first peak pressure and maximum mean pressures are predicted within an accuracy of respectively 29 and 12%. The rate of pressure rise dp/dt is typically underpredicted by 90%. The first and second eigen frequency of the pressure wave phenomena are predicted within 6%, while their corresponding amplitudes were predicted within 98%. Overall, it can be concluded that the current model predicts the considered experiments well.

6. Summary and conclusions

Fig. 24. The variation of (a) the axial flame front position, (b) the pressure with time and the variation of (c) the amplitude with frequency calculated using simulations for hydrogen–air–CO2 –He for four different diluent concentrations. The results corresponds to CFD simulations.

To summarize, the simulations very well reproduce the experimental trends of a decrease in the flame velocity, the intermediate peak pressure, the slope of the first dp/dt, and the maximum value of the mean pressure with an increase in the dilution. Furthermore, it is observed in the experiments that an increase in diluent mole fraction has almost no affect on the values of the first and second eigen frequency. This trend is also reproduced well in the simulations. This trend is consistent with the trends observed in Fig. 13. In addition, similar trends are observed in the simulations with H2 O as diluent as observed in the CO2 –He diluent. Furthermore, the flame dynamics predicted by H2 O and CO2 –He diluents are similar. However, the values of the maximum mean pressure obtained from H2 O simulations are somewhat lower than the values obtained from simulations corresponding to CO2 –He diluent. This trend is consistent with the trends presented in Section 4. 5.6. Quantitative comparison In order to further conclude about the predictive capability of our combustion model for the present experiments, it is necessary to make quantitative comparisons of the simulation and experimental results. For that purpose, we made a quantitative comparison of following parameters (a) intermediate peak pressure, (b) the maximum mean pressure, (c) the first rate of pressure rise (dp/dt), (d) the second rate of pressure rise (dp/dt), (e) the first eigen frequency, (f) the second eigen frequency, and (g) the amplitudes corresponding to these frequencies. These parameters have been selected because of their importance to the structural integrity of the containment. Table 2 shows a detailed comparison of these

The CFD based method described in our previous article (Sathiah et al., 2012b) is further extended to simulate turbulent flame propagation in hydrogen–air–diluent mixtures. The method consists of solving the density-averaged Navier–Stokes equations with the standard k −  turbulence model and an extended TFC combustion model. Validations of this TFC combustion model together with adaptive mesh refinement has been performed using three experiments performed in the ENACCEF facility. Namely, the runs RUN 153-10, 153-20 and 153-30 with hydrogen–air mixtures with three different CO2 –He dilution concentrations of 10, 20 and 30 vol.% have been used for the validation purposes. As a first step, the physical properties are compared for hydrogen–air–CO2 –He and hydrogen–air–H2 O mixture to assess whether CO2 –He can be used as a diluent instead of H2 O. In the next step, validation of the extended TFC combustion model was performed. For each validation case, a detailed mesh sensitivity study was performed by changing the number of levels of adaptive mesh refinement. For each diluent concentration, simulations have been performed for hydrogen–air–CO2 –He and hydrogen–air–H2 O mixtures. The main conclusions of the present work are: • The experiments show a decrease in the flame velocity, the intermediate peak pressure, the first rate of pressure rise dp/dt, and the maximum value of the mean pressure with an increase in the CO2 –He dilution. This decrease in the considered parameters is captured well in the simulations. In the experiments, an increase in the dilution has practically no effect on the first and the second eigen frequency. This trend is also reproduced well in the simulations. Furthermore, in the experiments, increase in the dilution from 0 vol.% to 20 vol.% dilution first decreases the amplitude of the residual pressure oscillations, while further increase in the dilution to 30 vol.% increases this amplitude significantly. In contrast, the CFD model predicts a systematic decrease in the amplitudes of the pressure fluctuations with an increase in dilution. This systematic decrease in the CFD analyses is consistent with the analysis of the physical properties presented in Section 4. • The simulation results of the flame dynamics, the intermediate peak pressure, and the slope of pressure with time for the different diluents (H2 O and CO2 –He) are very similar. However, the maximum mean pressure obtained by the simulation with CO2 –He as diluent is slightly higher in comparison to H2 O dilution. This is consistent with the trend observed from the analysis of the physical properties in Section 4. • The maximum value of the mean pressures and the intermediate peak pressures were predicted respectively within 12 and 29% accuracy. The rate of pressure rise dp/dt was typically underpredicted within 15–90%. The first and second eigen frequencies of the residual pressure wave phenomena were predicted within 6%. Therefore, it could be overall concluded that the presented CFD modeling predicts the considered ENACCEF experiments well.

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042

G Model NED-7927; No. of Pages 15

ARTICLE IN PRESS P. Sathiah et al. / Nuclear Engineering and Design xxx (2014) xxx–xxx

7. Future step In order to apply our combustion modeling to full scale containment applications, further extension of our combustion modeling to simulate non-uniform hydrogen–air mixture must be performed. Three experiments RUNs 765 and 736 as performed (see Table 1) in the ENACCEF facility (Chaumeix and Bentaib, 2011) can be used to validate the extended CFD modeling. The results of this validation will be presented in a follow-up paper. 8. Additional remarks The analyses presented in this paper also form part of a benchmark exercise. This benchmark exercise has been conducted with a number of participants within the SARNET-2 EU project. The results of the benchmark exercise will be presented in a complementary paper. Acknowledgment The work described in this paper is funded by the Dutch Ministry of Economic Affairs. References ANSYS FLUENT, 2008. Fluent 12.0. ANSYS-Fluent Inc., Lebanon, NH. Bentaib, A., Caroli, C., Chaumont, B., Chevalier-Jabet, K., 2010. Evaluation of the impact that pars have on the hydrogen risk in the reactor containment: methodology and application to PSA level 2. Science and Technology of Nuclear Installations. Bentaib, A., Chaumeix, N., 2012. SARNET H2 combustion benchmark diluent effect on flame propagation blind phase results. Tech. Rep., IRSN. Beyler, C.L., Custer, R.L., Walton, W.D., Watts, J.M.J., Drysdale, D., Hall, J.R.J., Dinenno, P.J., 2002. SFPE Handbook of Fire Protection Engineering. National Fire Protection Association. Blazek, J., 2005. Computational Fluid Dynamics: Principles and Applications. Elsevier Science Ltd. Bleyer, A., Taveau, J., Djebaïli-Chaumeix, N., Paillard, C., Bentaib, A., 2012. Comparison between FLACS explosion simulations and experiments conducted in a PWR Steam Generator casemate scale down with hydrogen gradients. Nuclear Engineering and Design 245, 189–196. CANTERA, 2011. Objected-oriented software for reacting flows. http://www.cantera. org/index.html/ Chaumeix, N., Bentaib, A., 2010. SARNET H2 combustion benchmark. Tech. Rep., Bureau de Physique des Accidents Graves, IRSN. Chaumeix, N., Bentaib, A., 2011. ISP 49-specification of ENACCEF test flame propagation in a hydrogen gradient. Tech. Rep., Bureau de Physique des Accidents Graves, IRSN. Fukushima, 2011. Fukushima Daiichi nuclear disaster. http://www.en. wikipedia.org/wiki/Fukushima$ $Daiichi$ $nucl%ear$ $disaster

15

GASEQ, 2011. A chemical equilibrium program for windows. http://www.c. morley.dsl.pipex.com/ Gubba, S.R., Ibrahim, S., Malalasekera, W., Masri, A., 2009. An assessment of large eddy simulations of premixed flames propagating past repeated obstacles. Combustion Theory and Modelling 13 (3), 513–540. Houkema, M., Komen, E., Siccama, N., Willemsen, S., 2003. CFD analyses of steam and hydrogen distribution in a nuclear power plant. In: The 10th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-10), Seoul, Korea, October 5–9. ISP-49, 2011. International standard problem ISP-49 on containment thermal hydraulics. Tech. Rep., Nuclear Energy Agency. Lamoureux, N., Djebaili-Chaumeix, N., Paillard, C., 2003. Laminar flame velocity determination for H2 –air–He–CO2 mixtures using the spherical bomb method. Experimental Thermal and Fluid Science 27 (4), 385–393. Liu, D., MacFarlane, R., 1983. Laminar burning velocities of hydrogen–air and hydrogen–air–steam flames. Combustion and Flame 49 (1–3), 59–71. Movahed-Shariat-Panahi, M., 2012. Recommendation for maximum allowable mesh size for plant combustion analyses with CFD codes. Nuclear Engineering and Design 253, 360–366. Poinsot, T., Veynante, D., 2001. Theoretical and Numerical Combustion. Roe, P., 1986. Characteristic based schemes for the Euler equations. Annual Review of Fluid Mechanics 18, 337–365. Roe, P.L., 1981. Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics 43, 357–372. Sathiah, P., Komen, E., Roekaerts, D., 2012a. The role of CFD combustion modeling in hydrogen safety management – I: validation based on small scale experiments. Nuclear Engineering and Design 248, 93–107. Sathiah, P., vanHaren, S., Komen, E., Roekaerts, D., 2012b. The role of CFD combustion modeling in hydrogen safety management – II: validation based on large scale experiments. Nuclear Engineering and Design 252, 289–302. Stempniewicz, M., 2009, April. SPECTRA Sophisticated Plant Evaluation Code for Thermal-Hydraulic Response Assessment Version 3.60, August 2009, Volume 1 – Program Description, Volume 2 – User’s Guide, Volume 3 – Subroutine Description, Volume 4 – Verification and Validation. NRG Report K5024/10.101640, Arnhem. Summers, R.M., Cole, R.K., Boucheron, E.A., Karmel, M.K., Dingman, S.E., Kelly, J.E., 1991. MELCOR 1.8.0: A Computer Code for Nuclear Reactor Severe Accident Source Term and Risk Assessment Analyses. Nuclear Regulatory Commission, Washington, DC. TECDOC, I., 2011. Mitigation of hydrogen hazards in severe accidents in nuclear power plants. Tech. Rep., IAEA-TECDOC-1661, IAEA. Visser, D., Agostinelli, G., Siccama, N., Komen, E., 2013. Hydrogen risk assessment – hydrogen distribution and mitigation. In: NURETH-15, Pisa, Italy, May 12–15. Visser, D., Houkema, M., Siccama, N., Komen, E., 2012a. Validation of a FLUENT CFD model for hydrogen distribution in a containment. Nuclear Engineering and Design 245, 161–171. Visser, D., Komen, E., Houkema, M., Siccama, N., Kyttälä, J., Huhtanen, R., Takasuo, E., 2009. FLUENT calculations of the hydrogen distribution in a containment during the OECD-NEA THAI HM-2 experiment. In: The 13th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-13), Kanazawa City, Ishikawa Prefecture, Japan, September 27–October 2, 2009. Visser, D., Siccama, N., Komen, E., 2012b. Hydrogen risk assessment – hydrogen distribution and mitigation analyses with CFD. In: TOPSAFE 2012 Conference, Helsinki, Finland. Zimont, V.L., 1979. Theory of turbulent combustion of a homogenous fuel mixture at high Reynolds number. Combustion Explosions and Shock Waves 15, 305–311.

Please cite this article in press as: Sathiah, P., et al., The role of CFD combustion modeling in hydrogen safety management – III: Validation based on homogeneous hydrogen–air–diluent experiments. Nucl. Eng. Des. (2014), http://dx.doi.org/10.1016/j.nucengdes.2014.05.042