Improvements in a CFD code for analysis of hydrogen behaviour within containments

Nuclear Engineering and Design 237 (2007) 627–647

Improvements in a CFD code for analysis of hydrogen behaviour within containments J.M. Mart´ın-Valdepe˜nas ∗ , M.A. Jim´enez, F. Mart´ın-Fuertes, J.A. Fern´andez Universidad Polit´ecnica de Madrid, College of Industrial Engineering, Jos´e Guti´errez Abascal 2, E-28006 Madrid, Spain Received 12 May 2006; received in revised form 6 June 2006; accepted 7 September 2006

Abstract The use of CFD codes for the analysis of the hydrogen behaviour within NPP containments during severe accidents has been increasing during last years. In this paper, the adaptation of a commercial multi-purpose code to this kind of problem is explained, i.e. by the implementation of models for several transport and physical phenomena like: steam condensation onto walls in presence of non-condensable gases, heat conduction, fog and rain formation, material properties and criteria for assessing the hydrogen combustion regime expected. The code has been validated against several experiments in order to verify its capacity to simulate the following phenomena: plumes, mixing, stratification and condensation. Moreover, two tests in an integral large enough experimental facility have been simulated, showing that the well-mixed and stratified conditions of the test were reproduced by the code. Finally, an example of a plant application demonstrates the ability of the code in this kind of problems. © 2006 Elsevier B.V. All rights reserved.

1. Introduction The systems and components of the Nuclear Plants were designed under safety criteria to avoid the maximum foreseen accidents. However, with very low probability, it has been demonstrated the occurrence of the so-called “severe accidents”, beyond the design basis. The TMI-2 accident (Broughton et al., 1989) showed the feasibility of this kind of accidents. Several research programs were started from this event for exploring the phenomena associated. Moreover, different studies, as the Probabilistic Safety Analysis (PSA), have been demanded by the regulatory bodies for evaluating and reducing the risk associated to these accidents non-considered at the design basis. These efforts have involved the main nuclear countries in the world. One of the main contributors to the containment early failure during a severe accident is associated to the presence of hydrogen within the containment building (Breitung, 1997). The hydrogen released from the reactor coolant system to the containment free volume, could mix or accumulate at different parts of the containment. If composition of the hydrogen–steam–air

∗

Corresponding author. Present address: C/ Pedro Justo Dorado Dellmans 11, E-28040 Madrid, Spain. Tel.: +34 91 330 0838. E-mail address: [email protected] (J.M. Mart´ın-Valdepe˜nas). 0029-5493/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2006.09.002

mixture lay within a certain limits, the combustion would occur. The combustion regime is a function of several variables as mixture composition, confinement and geometry, turbulence, etc. The pressure developed during the combustion process, which could threaten the containment integrity, depends on the regime and the amount of hydrogen burned. Therefore, the deep understanding of the hydrogen distribution throughout the containment at the accident time is very important. During the TMI-2 accident, 460 kg of hydrogen were released through the oxidation of 45% zirconium in the core. Most of the hydrogen was burned when its averaged concentration was 7.9 vol.%, and 35% of steam likely initiated by a spark in the relief tank room. The pressure rose at 3 bar maximum without any significant damage to the containment (Henrie and Postma, 1983, 1987). During last years, several studies and experimental research programs haven been performed around the world in the field of the hydrogen explosion risk in containment. The most relevant ones have been conducted by FZK (Fortschungzentrum Karlsruhe, Germany), KI (Kurchatov Institute Moscow, Russia), NUPEC (Nuclear Power Engineering Corporation, Japan), CEA (Commissariat a` l’Energie Atomique, France), NRC (Nuclear Regulatori Comision, USA), etc. (CSNI, 1999, 2000). These programs and several analyses performed have permitted the implementation of hydrogen risk mitigation measures at the current plants and their consideration in new designs.

628

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

Nomenclature Am,j B cp dAB dJ D Ea f f(Yg ) g Gr h hfg hfg H i iλ iσ k L m ˙ M n Nu p Pr PrT q R Ra Re t T um U V Vi x X Y z

coefficient j of the property m body force (kg/(m2 s2 )) specific heat (J/(kg K)) collision diameter (m) nozzle inner diameter (m) characteristic cloud size (m) activation energy (K) film thickness (m) degradation factor acceleration of gravity (m/s2 ) Grashof number (Gr = gρ∞ (ρi − ρ∞ )L3 )/μ2 heat transfer coefficient (W/(m2 K)) latent heat of vaporization (J/kg) latent heat of vaporization corrected by the subcooling effect (J/kg) total enthalpy (J/kg) static enthalpy (J/kg) DDT index FA index turbulence kinetic energy (m2 /s2 ) surface length (m) mass flux (kg/s) molecular weight (kg/kmol) maximum number of single gas species in the mixture Nusselt number pressure (Pa) Prandtl number (Pr = cp μ/k) turbulent Prandtl number heat flux (W/m2 ) universal gas constant (8.314 J/(mol K)) Rayleigh number (Ra = Gr × Pr) Reynolds number time (s) temperature (K) averaged velocity in the film (m/s) velocity vector (m/s) cloud volume (m3 ) volume of the node i normal to wall coordinate (m) molar fraction mass fraction axial coordinate (m)

Greek letters β thermal expansion (K−1 ) γ surface tension (N/m2 ) Γ diffusion coefficient (m2 /s) film thickness (m) δf ε turbulence dissipation rate (m2 /s3 ) θ wall inclination angle (rad) λ mixture detonation cell size (m) λK thermal conductivity (W/m) μ dynamic viscosity (Pa s)

μT ν ξ ρ σ σ* φ ΩAB

turbulent viscosity (Pa s) specific volume (m3 /kg) generic variable to be averaged density (kg/m3 ) expansion ratio critical expansion ratio equivalence ratio (φ = pH2 /2pO2 ) collision integral coefficient

Subscripts A specie A in the mixture b burned mixture B bulk specie (air) cd condensation cv convection f condensate liquid film fo fog g non-condensable gas H2 hydrogen H2,dry dry hydrogen H2 O steam i node of the room or cloud I interface l liquid phase m gas mixture Nu pure vapour condition O2 oxygen ra rain rd radiation ref reference sat saturation condition t total u unburned mixture w wall ␦ diffusion boundary layer ∞ bulk condition Superscripts = horizontal surfaces ave averaged DWN downward propagation INERT steam inert limit LFL lower flammability limit LOW lower condenser MED medium condenser T transposed UFL upper flammability limit UP upper condenser UPW upward propagation

Some examples are the implementation of passive autocatalytic recombiners (PARs) in PWR plants and inertisation and ignition devices in some BWRs. On the other hand, the hydrogen explosion risk is also a relevant safety topic in other fields. Some current industrial

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

applications as petrochemical, microchips and other make use of hydrogen as a primary constituent. Moreover, nuclear fusion reactors, in process of prototype design now, could present accidental sequences with important hydrogen releases (ITER, 2001). Finally, the future use of hydrogen as an energy carrier will increase the need of studying its safety management in different applications, in order to maintain the same risk levels than the fossil fuels which should be replaced (EUR, 2003). The analysis of hydrogen distribution and combustion in containment has to be faced in several kinds of approaches. Firstly, there are some phenomena of small time and spatial scale, such as buoyancy plumes, local accumulations or stratification, which need to be simulated in detail to find the required accuracy. For these phenomena, Computational Fluid Dynamics (CFD) codes are suitable. On the other hand, the long time of postulated accidental sequences, the large volume of the containments and the importance of the evolution during the long-term imply the use of more simplified codes, the Lumped Parameter category (LP). The Severe Accident Group at the Chair of Nuclear Technology (Polytechnic University of Madrid) has developed a methodology combining both kinds of codes for studying the hydrogen problem (Jim´enez et al., 2003a). This methodology consists of simulating a large enough number of sequences with a LP code to cover the different hydrogen release conditions expected during the accident and applying the CFD code to a selected number of sequences or short time windows in order to get an approach to the small-scale phenomena but avoiding an excessive computational cost. As the LP tool, this methodology applies the MELCOR code (Gauntt et al., 2000). Several institutions have developed specific CFD codes to cope with the hydrogen problem. The GASFLOW code (Travis et al., 1998), developed by Los Alamos and FZK, and the CAST3M code by CEA (Paill`ere et al., 2003), are currently used for deterministic analysis (Royl et al., 1996, 1997, 1999, 2000, 2001, 2002). The main feature of these codes is that the specific models for the hydrogen problem have been implemented during the development process. There is another group of CFD codes, the so-called multi-purpose, that have been developed by a commercial branch for several kinds of industrial applications. These codes usually feature a very large variety of models for different applications and a major development in the interface programs. Their main disadvantage is the lack of specific models for the analysis of hydrogen behaviour in containment. In this paper, the adaptation of a commercial CFD code, CFX4 (AEA, 2001), for the analysis of the hydrogen behaviour within containment is described. The main models introduced are: steam condensation onto walls in presence of non-condensable gases, heat conduction through the structures, fog and rain formation and specific models for evaluating material properties of the gases present in the containment atmosphere. Finally, this problem is not completely characterised just with the hydrogen distribution analysis within containment. The threat to the containment integrity is a consequence of the hydrogen combustion and the pressure set by the explosion regime developed. The deep analysis of the burning process is very complex, even for a CFD code. It needs specific and very detailed models to simulate, e.g. chemical reactions,

629

turbulence interaction with structures, etc. Moreover, the scale of this process implies very large both computational grids and simulation times: cell sizes below 0.01 m and time steps below 0.01 s (Breitung, 1997). Finally, some of the regimes could evolve to other ones, as by experiencing deflagrationto-detonation transition, which implies an increased level of complexity in model and analysis. Because of this complexity and the large uncertainties, the numerical simulation of the combustion processes has been excluded of our methodology. Instead, a prediction of the combustion regimes expected may be achieved by applying the widely known combustion criteria (CSNI, 2000), developed by FZK and KI, as functions of the mixture properties and geometry conditions. Code and implemented model descriptions are exposed in next chapters. Moreover, the validation tasks of the resulting code are shown in other chapter. These validations have been performed in a twofold strategy: separate-effect experiments and an integral, large enough, experimental facility. The paper is completed with an example of plant application, which demonstrates the capacity of the code to deal with the problem of hydrogen behaviour within the containment. 2. The CFD code The CFX-4 code (AEA, 2001) is a multi-purpose code developed for several industrial applications, which includes models for different kind of flows: single-phase, turbulence, multicomponent, multi-phase, combustion, etc. An important feature is the USER FORTRAN facility, which consists of a group of subroutines that can be accessed by the user in order to modify or to include own models. The code solves numerically the fluid dynamic transport equations by using the finite-volume approach in a structured multi-block mesh. Moreover, it incorporates turbulence closure models to solve the turbulence effects in reasonable computing times. The analysis of hydrogen behaviour within the containment, to be performed with this code, considers the following flow features: single-phase, compressible, multi-component, turbulent and non-reactive. The set of equations to be solved with CFX-4 are the Navier–Stokes equations in their conservation form: The continuity equation: ∂ρ  =0 + ∇ · (ρU) ∂t

(1)

The momentum equation: ∂ρU  ⊗ U)  = (ρ − ρ0 )g + ∇ · (μ + μT )((∇ · U)  + ∇ · (ρU ∂t   2 T  + (∇ · U) ) + ∇ · p + ρk (2) 3 The energy equation:   ∂ρH μT ∂p  ∇i = + ∇ · (ρUH) − ∇ · λK ∇T + ∂t PrT,H ∂t

(3)

630

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

The total enthalpy is expressed in terms of the static enthalpy and kinetic energy as: 1 H = i + U2 2

(4)

The transport equations of the variable mass fractions (YA ) are:     ∂ρYA μT  A − ΓA + ∇YA = 0 (5) + ∇ · ρUY ∂t PrT,A The code has several turbulence models, but in this analysis the standard k–ε model has been used. This model assumes that turbulence viscosity is linked to turbulence kinetic energy k and turbulence dissipation rate ε via a relation (μT = C␮ ρ(k2 /ε)). The values of k and ε come directly from their differential transport equations:    ∂ρk μT  −∇ · ∇k =PT + GT − ρε + ∇ · (ρUk) μ+ ∂t PrT,k (6)

∂ρε  −∇ · + ∇ · (ρUε) ∂t



μT μ+ PrT,␧



 ∇ε

ε ε2 = C1 (PT + C3 max(GT , 0)) − C2 ρ k k

(7)

PT is the shear production term and GT is the buoyancy production term: T

  + (∇ · U)  ) PT = (μ + μT )∇ · U((∇ · U) 2   + ρk) − ∇ · U((μ + μT )∇ · U 3 GT = −

μT g · ∇ρ PrT,␳

The model constants are given by: C␮ = 0.09, C1 = 1.44, C2 = 1.92 and C3 = 1. The turbulence Prandtl numbers are: PrT,k = 1, PrT,␧ = 1.217 and PrT,H = PrT,A = PrT,␳ = 0.9. These expressions are completed by adding two algebraic equations—the equation of state system and the constitutive equation: ρA =

pMA , RT

0≤A≤n

(8)

n  YA 1 = ρ ρA

(9)

A=0

h=

n 

 YA hA =

A=0



−

1 Tref

 0

n 

A=0 Tref

1 YA T



T

 

cpA (T ) dT



0

cpB (T  ) dT 

 (10)

being the last term of Eq. (10), the reference enthalpy, i.e. enthalpy of the bulk fluid at the reference temperature (Tref ) where static enthalpy is defined to be zero. Several discretisation schemes for the advection terms are available in the code. The simulations presented in this paper have been performed with a “Hybrid” scheme that combines an upwind and a central-difference scheme as a function of the Peclet number in the mesh (AEA, 2001). As general rules for meshing, in relation with the condensation models implemented in the code, the treatment of the boundary layer is of particular interest. The default wall function for the k–ε model has been used in this study. For this, the first grid point (adjacent to the wall) must lie within the logarithmic part of the boundary layer. A profile including the viscous sub-layer and the logarithmic region is applied therein. In case of a finer grid near the wall, this restriction is violated and the wall function accuracy will deteriorate. 3. Additional model implementation Steam and hydrogen, generated in the reactor coolant system during the in-vessel phase of a severe accident, are released to the containment atmosphere through the break point. The formation of expanded light clouds of released gases establish jet and buoyant plume structures as a result of inertial and buoyant forces, which make light gas rise up. When the plume impinges the uppermost containment, it spreads throughout the dome ceiling and stagnates there. Depending on the release location and the geometrical aspect ratio (slenderness) of the building, the inertia forces would be able to drive the atmosphere towards either wellmixed or stratified conditions. In the medium- and long-term, other mixing phenomena could arise and change the atmosphere conditions, mainly heat and mass transfer processes (convection and condensation). Moreover, containment engineered systems could change these conditions, normally by enhancing the mixing processes, i.e. fan coolers, internal/external sprays, rupture disks and passive autocatalytic recombiners. All of these phenomena govern the hydrogen distribution status within the containment: well-mixed, stratified, locally gathered, sedimented, etc. (CSNI, 1999). Several containment phenomena can be simulated by the code. Some examples are summarised in Section 4 below. However, there are some specific phenomena which have not been included by the code developers and are very relevant for the hydrogen problem. The most important phenomenon is the condensation onto walls. It controls the steam presence in the containment in the medium- and long-term, leading to inert or flammable conditions. Moreover, some of the mixing/gathering processes are driven mainly by condensation. The presence of non-condensable gases in the containment degrades the steam condensation on the walls. In order to permit adequate simulation of the fluid-dynamics in the containment, other models have been implemented as: heat conduction through the structures, fog and rain formation and specific models for estimating the material properties of the atmosphere gases. Finally, the assessment of the combustion regime is essential information for these studies, hence the combustion criteria have been implemented in

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

the code as well. Next, the formulation of these models and some information of the implementation procedure are summarised, details may be found in Mart´ın-Valdepe˜nas (2004).

Different condensation models have been implemented and linked to the CFX code. A comparison of these models and their performance against several experimental databases are explained in Mart´ın-Valdepe˜nas et al. (2005). In the current paper, only the most suitable model for the subsequent application is detailed. The first assumption is the steady-state behaviour at the interface: the heat power conducted through the liquid film attached to the wall is equal to that coming from the atmosphere side by convection, radiation and condensation: qt = qIw = q∞I = qcv + qrd + qcd

(11)

The radiation heat transfer can be neglected since the bulk temperature is relatively low during most of the anticipated sequences (Terasaka and Makita, 1997), being not considered the combustion processes in these studies. The previous equation can be expressed in terms of the heat transfer coefficients and temperature differences as follows: qt = ht (T∞ − Tw ) = hf (TI − Tw ) = hcv (T∞ − TI ) + hcd (T∞,sat − TI )

(12)

The models described below allow calculating these heat transfer coefficients, apart from the convective heat which is directly calculated by the code by using wall functions. 3.1.1. Film models The heat transfer coefficient through the condensate liquid film depends on: (1) film flow regime, (2) surface orientation and (3) temperature profile through the liquid film. Considering a leaning wall, the film heat transfer coefficient for laminar flow is calculated by Nusselt theory (Collier, 1972): hf =

ρl (ρl − ρ∞ )g sinθ hfg (λK,l )3 4μl z(TI − Tw )

expressed in terms of the Reynolds number for the film: Ref =

4ρl um δf μl

(15)

The Nusselt number in the film may be defined as follows:

3.1. Condensation models



631

1/4

Nuf =

(3/4)hf [μ2l /ρl (ρl − ρ∞ )g sin θ] λK,l

1/3

(16)

Using the expressions of Ref and Nuf , Eq. (13) may be rearranged (Collier, 1972): −1/3

Nuf = 1.47 Ref

(17)

this equation may be used in laminar wave-free flows (Ref ≤ 30). In laminar wavy films (30 ≤ Ref ≤ 1800) Kutateladze (Incropera and DeWitt, 1996) recommended: Nuf =

Ref 1.08 Re1.22 − 5.2 f

(18)

and, for turbulent films (Ref ≥ 1800), Labuntsov (Incropera and DeWitt, 1996) suggested: Nuf =

Ref −0.5 8750 + 58 Prl (Re1.22 f

− 253)

(19)

The previous models are valid for vertical or inclined walls. When the surface is nearly horizontal other effects such as gravity and capillarity become important and these models may not be valid. In these situations, specific models for horizontal walls have been implemented. The Nusselt number is modified to take into account the surface tension of the liquid: Nu= f =

hf [γl /(ρl − ρ)g cos θ]1/2 λK,l

(20)

and, for the Rayleigh number for the film: Raf= =

g cos θ(ρl − ρ∞ )hfg L3 μl λK,l (TI − Tw )

(21)

For downward-facing surfaces and laminar flows (Raf= < 10 ), the correlation is (Gauntt et al., 2000): 8

(13)

The liquid properties are evaluated at the intermediate temperature between wall and interface (Tf = (1/2)(Tw + TI )), hfg is evaluated at TI and ρ∞ at bulk conditions. Notice that Eq. (13) considers a static liquid film which drains along all the nodes adjacent to walls. The condensate subcooling and temperature jump across the film are considered by:

 cpl (TI − Tw ) (14) hfg = hfg 1 + 0.68 hfg The liquid film at the wall surface may drip toward the lower part by gravitational forces, and the film may develop a rippled structure or turbulent flow. The transition criterion may be

0.2

6 = Nu= f = 0.6[max(10 , Raf )]

(22)

and, for turbulent film regimes (Raf= ≥ 108 ): 0.19

10 = Nu= f = 0.72[max(10 , Raf )]

(23)

For upward-facing surfaces, the Clifton–Chapman model (Chapman, 1984) considers the motion only by the hydrostatic force:  1/4   gρl2 hfg L3 λK,l (TI − Tw ) = Nuf = 2.43 (24) F μl hfg μl λK,l (TI − Tw ) being the function F approximately equal to 0.3 in the containment conditions.

632

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

Table 1 Experimental conditions, test results and comparison with the condensation model results at the Dehbi’s tests Test

T∞ (◦ C)

Tw (◦ C)

p (bar)

YAIR

hEXP (W/(m2 K)) t

hCFD (Eq. (26)) (W/(m2 K)) t

A38 A30 A25 B39 B33 B28 C25 C16 C09

101 94 79 125 113 85 137 127 95

72 65 60 90 88 62 100 82 60

1.5 1.5 1.5 3 3 3 4.5 4.5 4.5

0.33 0.56 0.80 0.3 0.59 0.85 0.35 0.58 0.88

788.4 461.2 189.7 981.8 502.1 141.3 1171.5 606.2 178.5

616.3 288.7 89.7 584.8 271.4 86.6 587.8 297.6 64.3

3.1.2. Heat and mass transfer with non-condensable gases The degradation effect of the non-condensable gases accumulating near the condensate film was considered by different models. Three models have been introduced based on: (1) an experimental correlation, (2) the heat and mass transfer analogy, (3) the diffusion boundary layer and (4) a variation of the third model, by replacing the convective heat transfer coefficient calculated by the CFD code with a heat transfer correlation. The first model works well in several experimental situations, with rather coarse grids and spends one-half of the CPU time the others do. The second one uses the Chilton–Colburn correlation (Collier, 1972). This model has shown the worst behaviour in most conditions, especially at high atmosphere pressure and with condensation onto horizontal walls. Moreover, it expends more computing time and needs finer meshes than the previous one. The last two models have shown similar behaviour when the CFD code estimated the convective heat transfer process accurately. However, in the case of horizontal surfaces, the wall law works badly and the condensation process is sub-estimated. The sophistication of these models imposes larger costs in computing time. Based on its good performance and low computational cost, the first model is usually applied by the authors in plant applications and is the only one presented here, more details can be found in Mart´ın-Valdepe˜nas et al. (2005). The correlation-based model was developed by Terasaka and Makita (1997). The steam condensation in a pure vapour atmosphere is calculated with the Nusselt theory, using Eq. (13), but under saturation conditions: 

qNu

ρl (ρl − ρ∞ )g sin θ hfg (λK,l )3 = 4μl z(T∞,sat − TI )

1/4 (T∞,sat − TI )

(25)

The condensation heat flux with non-condensable gases is calculated with the previous equation but affected by a factor reducing its value, the “degradation factor”: qcd = f (Yg )qNu

(26)

The factor was fitted by Terasaka and Makita to different experimental databases: Asano et al. (1979), PHEBUS (1997) and 2D numerical results by Kim and Corradini (1990) in the mass fraction range from 0 to 1. The resulting expression

was: f (Yg ) = [1 − 0.964Yg + 4.989Yg2 − 4.135Yg3 ] ·

1 − Yg 1 + 15.48Yg (27)

Finally, the steam mass flow towards the liquid film is calculated from: qcd m ˙ cd =  (28) hfg The correlation Eq. (27) was compared by Terasaka and Makita (1997) with the UCB database (Vierow and Schrock, 1991), giving good results in the whole range of data. The quite different experimental conditions used to fit the correlation provide the model with a wide applicability range. However, in this experiment, some scattering of data was identified and not reproduced by the correlation. The reason of these lacks of the models was not explained in Terasaka and Makita (1997). However, this effect has been pointed out in the validation exercise, presented below, as a pressure effect and it has been explored through the following factor (Mart´ın-Valdepe˜nas et al., 2005): 2/3  qcd p = f (Yg ) (29) qNu 1.15 × 105 the results of comparison with the experimental data are depicted below in Table 1. This factor reflects the pressure influence affecting the model, as was indicated by the mechanistic analysis of Herranz et al. (1998): ht ∼ Gr1/3 ∼ (ρg,∞ )2/3 ∼ p2/3 . 3.1.3. Iterative algorithm The simulation of the condensation process with these models needs to evaluate the properties and the equations at the interface temperature (TI ), which is an unknown of the problem. Therefore, an iterative procedure is necessary to find the heat transfer coefficient. The algorithm and the models have been implemented in CFX-4 by using the so-called USRSRC userroutine, which permits to add source terms in the governing equations. The interface temperature is estimated in the first step of the algorithm. The properties and heat transfer coefficients of the film are calculated with this datum in each node adjacent to convective boundary walls. The convective heat transfer is

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

calculated by the code using the wall law. The condensation heat flux is then calculated, Eq. (26) or (29). Finally, a heat balance in the interface is performed, Eq. (12), and a new value of TI is found by: TI =

Tw hf + T∞ (hcv + hcd ) hf + (hcv + hcd )

(30)

If TI is close enough to the current temperature the algorithm is stopped, otherwise, the procedure is restarted with the new TI . The algorithm is run with a precision of 0.001 K. When the iterative procedure is converged mass, energy and steam sources/sinks are imposed onto the continuity, energy and steam equations in the node selected. This procedure is repeated in every node adjacent to the convective boundary walls and in each iteration at time step. 3.2. Additional thermalhydraulic models Other important thermalhydraulic phenomena nonreproducible by the code have been implemented via simplified models. They are: (1) heat conduction through the structures, (2) fog and rain formation and (3) specific models for estimating the material and mixture properties of the gases. Other models as sprays performance, pool formation, PARs actuation, etc., have not been implemented yet but could be incorporated in future versions of the code. 3.2.1. 1D heat conduction The main phenomena related to the long-term mixing process in containment are the heat and mass transfer within the containment walls, i.e. convection and condensation. This mechanism depends on the wall temperatures, based on heat conduction and storage in the structure. Although it is a 3D phenomenon, the use of a 1D formulation gives precision enough and simplifies the numerical array and the geometry construction. Moreover, it is known that the wall temperature is affected by the heat flux by condensation and convection more than by the transverse heat conduction. The transient 1D heat conduction equation in Cartesian co-ordinates is: ∂ρH ∂2 − 2 · (λK ∇T ) = 0 ∂t ∂x

(31)

This equation is discretised through either an implicit or an explicit scheme, as selected by the user. The USRBCS userroutine, which permits to specify boundary conditions at walls in each time step or iteration, has been used in this case. Therefore, the wall temperature used in each time step to calculate the convection and condensation heat fluxes is evaluated in each boundary node by the solution of the 1D heat conduction equation in each conducting structure node. The solution of Eq. (31) at the structure bound considers different kinds of boundary conditions: (1) a heat flux between fluid and wall by condensation and conduction, (2) a constant temperature and (3) an adiabatic condition. The user can select the best condition in each node.

633

3.2.2. Fog and rain formation The fog formation in closed systems happens when the vapour pressure is higher than the saturation pressure in the atmosphere. In such conditions, the steam condenses, forming little droplets without any cold surface. The steam will condense on these droplets and grow them up. There are several models to simulate these phenomena; the majority are based on the droplet kinetic theory via the Mason equation (Collier, 1972). These models are complex but need some experimental constants fitted to the experiments. In order to avoid this problem, a very simple model, used by other containment CFD code (Travis et al., 1998), has been selected: m ˙ fo = Cfo V (ρH2 O,sat − ρH2 O,∞ ) ρH2 O,sat (T∞ , psat (T∞ )) =

psat (T∞ )MH2 O RT∞

(32) (33)

being Cfo an experimental constant, which values usually 0.1 s−1 (Travis et al., 1998). The fog droplets are small enough to stand airborne as aerosols. They move in the containment bulk as other gases. When the droplet density is high, they coalesce to larger droplets, which cannot be airborne anymore and fall by gravity as rain if larger than a critical size. The model implemented makes use of a relaxation function similar to Eq. (32) and was taken from the same reference (Travis et al., 1998): m ˙ ra = Cra V (2 × 10−3 − ρra,∞ )

(34)

Being 2 × 10−3 kg/m3 the maximum value of fog density within the atmosphere and Cra the experimental constant, usually 0.01 s−1 (Travis et al., 1998). 3.2.3. Material and mixture properties The code CFX-4 permits, in its command file, to introduce constant properties (cp , μ, λK , Γ ) for each species. However, in a severe accident, the gas phase conditions could experience important changes: temperatures from 300 to 1000 K and pressures from 1 to 5 bar. Therefore, these variations have to be considered in the calculations. Moreover, some of the implemented models need these properties for each species or mixture. Some additional modules have been implemented with this aim. Properties of the single gas species have been calculated from the polynomial functions of temperature fitted by the NASA, compiled in the CEA-Program (Gordon and McBride, 1994). These functions have been added to the CFD code for the following gases: air, nitrogen, oxygen, hydrogen and helium and m represents μ and λK : ln m = Am,0 ln T + Am,1 T −1 + Am,2 T −2 + Am,3

(35)

cp = Ac,0 T −2 + Ac,1 T −1 + Ac,2 + Ac,3 T + Ac,4 T 2 R + Ac,5 T 3 + Ac,6 T 4

(36)

Properties of water (liquid and vapour) have been taken from the TABAGUA code developed by Centro At´omico de Bariloche

634

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

(CAB, 1988). The program uses functions to calculate cp , μ, λK , Tsat , Psat , σ l and β. The heat of vaporisation is taken form DIPPR database (AEA, 2001). The binary diffusion coefficient of the species A and B is calculated with the correlation of Wilke and Lee (Bird et al., 1960): T 3 ((0.001/MA ) + (0.001/MB )) −2 ΓAB = 1.88292 × 10 2 Ω pdAB AB (37) The models implemented in the CFD code for calculating the condensation process are based on the mixture properties, generally near the condensing wall. The mixture models used here have been taken from (Reid et al., 1988) and are calculated at T␦ . 3.3. Hydrogen combustion criteria During a severe accident, different hydrogen combustion regimes could happen, because the hydrogen concentration could reach until 20 vol.% and the steam concentration could range between 20 and 70 vol.% (CSNI, 2000). Moreover, the mixing and distribution processes could drive to either well-mixed or inhomogeneous conditions (stratification, local accumulation, etc.). If an ignition source were present the hydrogen could burn in different conditions, depending on the hydrogen and oxygen concentration and the flow geometry. The combustion regimes, from the point of view of the nuclear safety, are classified in CSNI (2000) as—(1) slow deflagration (SD): the flame velocity is lower than the reactant sound velocity (>200 m/s) and the pressure increase is of the order of the initial one. (2) Flame acceleration (FA): the flame velocity is higher than the reactant sound velocity, but lower than product sound velocity (>500 m/s) and the pressure increases some 10 times the initial one. (3) Deflagration-to-detonation transition (DDT): the flame velocity is higher than the product sound velocity (∼1200 m/s) and the pressure increases about thirty times the initial one. The simulation of the combustion process with a CFD code is a very complex task as was discussed at the end of the first section. Therefore, we have used three criteria for identifying the combustion regimes (CSNI, 2000), avoiding combustion simulations: (1) the flammability of the mixture (flammability criterion), (2) the possibility of FA (σ criterion) and (3) the possibility of DDT (λ criterion). 3.3.1. Flammability criterion The burning process occurs when there is an ignition source and the mixture composition lies within the flammability limits, i.e. the hydrogen/air/steam proportion is adequate for sustaining combustion. This study neglects the ignition source, because of the random feature of this process during an accidental sequence, and only focuses on the gas concentration. The factors that affect the flammability limits are: (1) hydrogen/air/steam concentration, (2) lean or rich hydrogen mixtures, (3) flame propagation orientation and (4) initial temperature conditions.

Fig. 1. Shapiro–Moffette ternary diagram showing the flammability limits for hydrogen/air/steam/mixtures at 300 K and 1 bar.

The flammability limits are not themselves a property of the mixture, but are very useful for engineering and safety studies. They are indicative of the mixture capacity to maintain the flame propagation after a high enough energy release by an ignition source (Stamps and Berman, 1991). Several experimental programs had evaluated these limits (Kumar, 1985; Marshal, 1986) for several temperatures. They usually are expressed in charts, as the widely known Shapiro–Moffette diagram (Fig. 1). Byun et al. (2000) proposed linear functions of the molar fraction of hydrogen, steam and temperature for expressing the flammability limits. From the data compiled by Stamps and Berman (1991), the following expressions are proposed: LFL-UPW XH ≥ 0.037 + 0.011XH2 O − 4.16 × 10−5 (T − 373) 2

(38) LFL-DWN XH ≥ 0.086 + 0.008XH2 O − 1.02 × 10−4 (T − 373) 2

(39) UFL XH ≤ 0.772 + 1.087XH2 O + 2.71 × 10−4 (T − 373) 2

(40)

INERT ≥ 0.63 + 3 × 10−4 (T − 373) XH 2O

(41)

3.3.2. σ criterion Once the mixture flammability is assessed with the previous criterion, the combustion regime is relevant information to evaluate the threat to the containment building during the accident by hydrogen combustion. The first step is checking the possibility of FA, by using the σ criterion developed by FZK and KI (CSNI, 2000). They performed several experimental series (tubes from 0.08 to 2.25 m in diameter, several blockage ratios and mixture concentrations) and reviewed the experimental data of the literature. They concluded that the FA does not depend on the geometrical scale or blockage conditions, but depends on the initial properties of the mixture (Breitung and Royl, 2000), if

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

enough turbulence is developed. The best parameter correlating the FA occurrence is the expansion ratio, called σ parameter: vb ρu σ= = ρb vu

(42)

The burned mixture is assumed through the adiabatic isobaric complete combustion of ideal gases. Dorofeev et al. (2001) have demonstrated that there is a critical value (σ * ) separating the mixtures with just SD apart from the mixtures undergoing FA. This value has been fitted by the following polynomial function (being φ between 0.1 and 7.5):  3  2   Ea Ea ∗ −5 Ea σ = 0.9 × 10 − 0.0019 + 0.1807 Tu Tu Tu + 0.2314

(43)

635

has the following size (Breitung and Royl, 2000): D = V 1/3

(47)

The λ parameter is a property of the mixture depending on pressure, temperature and gas concentrations. A wide experimental database exists for this property, hence KI have fitted different analytical functions by minimum squares (Dorofeev et al., 1996). The following, fitted under the uncertainty limits of the experiments, has been implemented:  log10 (λ) =

−1.13331 − 2.38970  +



4.59807 × 101 0.997468

((XH2 ,dry + 4.18215 × 102 )/T )

2

Dorofeev et al. (2001) proposed the following expression for the activation energy, which was fitted by minimum squares to experimental and thermodynamic data: Ea,H2 = 7.73 × 103 − 4.06 × 102 φ + 8.96 × 101 φ2 − 4.32 × 10−1 φ4 The FA index is: σ iσ = ∗ σ

(44)

(45)

When iσ ≥ 1.0, FA is possible and, if iσ < 1.0, FA is not possible. However, uncertainties in the model (8%) recommend reducing this threshold value to 0.92. 3.3.3. λ criterion DDT processes occur when the local core, due to the pressure wave, is amplified and propagated from the reaction zone to the surroundings of the non-perturbed mixture. A stable coupling between the shock wave and the reaction zone has to happen. From numerical simulations, Dorofeev et al. (1989) suggested that this process needs a minimum size of the reactive mixture. The subsequent experimental program performed and the reviews of the available data in the literature have confirmed this result (Breitung and Royl, 2000). The analysis showed (on facilities of different characteristic sizes) that the minimum mixture size for DDT must be above a certain value measured in units of the detonation cell size, known as the λ parameter. Hence, the minimum characteristic cloud size for DDT occurrence is 7λ (Dorofeev et al., 1996). The DDT index is therefore defined as: iλ =

D 7λ

(46)

When iλ ≥ 1, DDT is possible and, if iλ < 1, DDT is not possible. However, uncertainties in the model (43%) recommend reducing this threshold value to 0.57. The characteristic mixture size is calculated from flammable cloud sizes in rooms or mixtures of different sensitivity, which

+ 8.74995 × 10−4 (XH2 ,dry + 2.66646 × 10−2 T )  2

− 4.07641 × 10−2 (XH2 ,dry +2.66646×10−2 T )

× (1 + 4.65429 × 10−2 XH2 O + 3.59620 2 −7 2 3.31162 × 10 ×10 T (XH2 O ) ) T

1 × (p + 1.57650 × 10−1 ) 0.1 + 1.57650 × 10−1 − 8.42378(p − 0.1) + 2.38970

(48)

3.3.4. Implementation in the CFD code Criteria defined above depend on mixture properties (gas concentration, pressure and temperature) before the combustion happens, avoiding the full simulation of the reactive flow in safety analysis. Moreover, other geometrical properties are used in this criteria, they are: the limits of the room where the combustion could happen and the flammable cloud size. The three criteria have been implemented via the user-subroutine USRTRN, which post-processes the results of these variables at the nodes or group of nodes in each time step, performs a check and indicates whether the criteria are fulfilled. The equations (Eqs. (38)–(41)) are checked at each node of the fluid field. If the mixture lied within the flammability limits, the combustion would occur when an ignition source existed. The hydrogen combustion is discarded otherwise. Next step consists of checking the σ criterion, which should be evaluated in a group of nodes. This criterion, as implemented in the CFX code, is checked in two levels, called here as GLOBAL and LOCAL. The first level considers all the nodes in a containment room and applies Eqs. (42)–(45) to the averaged values of each variable in each room (CSNI, 2000). This parameter permits to evaluate effects of the non-flammable gas on the flammable gas in rooms where the mixture is unhomogeneous. Averaging is made in a volume basis for all the variables in the σ criterion (ξ) by using

636

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

the following expression:

NODES ξi V i ¯ξ = i NODES Vi i

(49)

The LOCAL level applies the criterion only to the flammable cloud, i.e. the variable average Eq. (49) is now applied only to the room nodes within the flammable limits, following Eqs. (38)–(41) of the flammability criterion. Even though this is not a rigorous application of the σ criterion, it is an indicator of the mixture sensitivity and a conservative way to reduce the uncertainties of the room limits and gradients. The λ criterion is also implemented in the CFD code in the two levels GLOBAL and LOCAL. The first one, in this case, considers all the nodes of the room which are within the flammability limits Eqs. (38)–(41). Variables in the selected nodes are averaged with Eq. (49) and the criterion is applied following Eqs. (46)–(48). The LOCAL level is applied to four clouds of different mixture sensitivities in each room, using the sane scheme as the GLOBAL one. The following clouds are averaged with Eq. (49): Cloud ‘‘tur’’ : 100XH2 + 100XH2 +

2.5 100XH2 O ≥ 6.5 58

Cloud ‘‘ace’’ : 100XH2 +

(50) 4. Validation

1 100XH2 O ≥ 9 and 26

1 1 100XH2 + 100XH2 O ≤ 1 70 62 Cloud ‘‘ddt’’ : 100XH2 +

enough separated from one another and candidates to accumulate important quantities of hydrogen, as dome, loop rooms, pressuriser room, annular regions, etc.

and

63.5 100XH2 O ≤ 72.5 58

(51)

1 100XH2 O ≥ 11 and 45

1 1 100XH2 + 100XH2 O ≤ 1 64 55 Cloud ‘‘det’’ : 100XH2 ≥ 15, 1 1 100XH2 + 100XH2 O ≤ 1 64 55

(52)

100XH2 O ≤ 30

Fig. 2. Limits of the clouds “inf”, “ace” and “det” at the Shapiro–Moffette diagram at 373 K.

and (53)

These clouds have been selected approximately from typical regions in the Shapiro–Moffette diagram at 373 K which would burn as turbulent flame, FA, DDT or on-set detonation, respectively. These regions are depicted approximately in Fig. 2. These criteria are applied to all the containment fluid nodes, but lumped in each room. These imply a previous selection of these rooms during the geometry construction. The CFX code has a feature, which permits to group some nodes, defining the so-called USER3D regions. Once the geometry is built up and previous to the mesh generation, these regions are selected with the pre-process program, and all the nodes in these regions are easily managed by the user subroutines. Criteria to select these regions consist on taking the main rooms of the containment,

This section is dedicated to the validation tasks of the developed code improvements. The aim is to evaluate and calibrate the code abilities to simulate the processes involved in the hydrogen mixing and distribution within large enclosures, as NPP containments. Two kinds of experiments have been simulated. In one hand, simplified experiments in small or medium scale facilities involve only one or a limited number of phenomena. On the other hand, integral experiments as a large-scale facility with rooms and compartments similar to actual PWR containments. Therefore, the most important mixing and distribution phenomena in the containment have been accounted for. In this section, the most representative examples in the extended validation exercise are summarised. Only the results of the optimal mesh for each simulation are presented, i.e. the ones balancing best between accuracy and computational cost. At these meshes, the numerical diffusion is reduced below the order of turbulent viscosity in order to avoid artificial mixing. Further details and additional benchmark exercises are introduced in Mart´ınValdepe˜nas (2004) and other references indicated below. 4.1. Basic experiments Several separate-effect experiments have been performed in order to check the behaviour and predictive abilities of the CFD code in simulating a jet, a light gas plume, natural convection heat transfer, stable/unstable semi-confined atmospheres, etc. (Mart´ın-Valdepe˜nas, 2004). The code demonstrated its capacity in simulating these separated phenomena with accuracy in most of the cases. Some medium-scale and more complicated experiments are presented herein; they have been selected to

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

Fig. 3. Comparison between calculated helium molar fractions and experimental results, in a vertical intermediate line and different times, at the helium experiment in a vented cavity.

demonstrate the code capacity in simulating the stratification phenomena coupled to the validation of the implemented condensation model. 4.1.1. Helium distribution in a vented cavity The experiment (Renson and Goldstein, 1987) consisted on a helium constant mass-rate release (4.6 × 10−3 kg/s) through a dJ = 0.1 m diameter nozzle, at the bottom of a cylindrical vessel 7 m3 (2.29 m high, 2 m diameter) volume. A vent opening at the top, kept the tank depressurised during the release phase (120 s). The calculation was performed with a 1200-cell 2D mesh and refining of the mesh at the lowermost cavity was found the optimal approach to the experiment. The code reproduces the turbulent circulation motion of the jet and the stratification of helium, which was accumulated in the uppermost half of the cavity. The experimental results, i.e. the stratification and the concentration gradients, were predicted accurately. The deviation on the helium concentration calculated, as is depicted in Fig. 3, was below 10% showing a delayed time evolution. 4.1.2. Hydrogen release and mixing in an enclosure Shebeko et al. (1988) performed a hydrogen distribution experiment for a subsonic release of hydrogen in a closed vessel. Hydrogen was released vertically upwards (3.7 × 10−4 kg/s) through a dJ = 0.01 m diameter nozzle at 1.4 m under the top of the vessel of 20 m3 (5.5 m high, 2.2 m diameter). The release period was 60 s, the injection was stopped afterwards and the hydrogen concentration was measured during 250 min. A 2D cylindrical mesh with 22,000 cells was considered the optimal one in this simulation. The mesh was refined at the upper part, near the release point. The code reproduced qualitatively the behaviour of the hydrogen at the experiment: initial stratification over the release point and subsequent homogenisation. However, the results, in comparison with the experimental data, were worse than in the previous benchmark (Fig. 4). At the vessel top the code overestimated the hydrogen in approximately 70% during all the experiment. On the other hand, the hydrogen below the release point was underestimated at the later phase of the test,

637

Fig. 4. Comparison between calculated hydrogen volume fractions and experimental results, in the vessel symmetrical axe and different times after the release phase, at the hydrogen release and mixing experiment in an enclosure.

scarcely a few hydrogen reach this area in the simulation. This deviation was based on the boundary condition of the simulation: adiabatic walls, which imposed a mixing process only by diffusion at the lower vessel. Being ΓH2 = 0.74 × 10−4 m2 /s, the mixing time would be around 7000 min, which is higher than the 250 min of the tests. This implies that other mixing phenomena, probably convection, participated in this experiment, but was not reported by the experimentalists (Gallego et al., 2005). 4.1.3. Condensation with non-condensable gases Several tests have been simulated in order to validate the new implemented condensation model in presence of noncondensable gases. Two series will be illustrated below. The first experimental series are that conducted by Dehbi et al. (1991) within a cylindrical vessel 5 m long and 0.45 m diameter. A condensing cylinder (3.5 m long and 0.0038 m in diameter) was placed inside the vessel, which was fully insulated so that condensation took place only onto the condenser surface. The steam was generated at the bottom of the vessel at a rate keeping the vessel pressure at a constant predetermined value. The steady-state total heat transfer coefficient (ht ) was measured in the test and has been used for comparison with simulation results. A 2D cylindrical mesh of 1500 nodes was considered as the optimal one for this experiment and the condensation model. Finer grids near the walls reproduced worse the condensation process because of the model dependence on the bulk properties. Nine tests were simulated, as shown in Table 1. The condensation model Eq. (26) sub-estimated the total heat transfer between 22 and 64%. The worst results were found at the higher pressure, lower non-condensable mass fraction and higher bulk temperature tests. From these results, the model is clearly independent on pressure. The main reason for that were the experimental databases used to fit this correlation (Uchida and PHEBUS), all of them performed at 1 bar of air. Therefore, the pressure influence was not considered in the model. This issue was already pointed out by Peterson (1996) for the Uchida model. Using the model modified with the pressure effect, Eq. (29), the deviations are clearly reduced, being the maximum one in the order of the experimental uncertainty (±25%) (Fig. 5).

638

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647 Table 3 Comparison of global variables with the condensation model results at the MICOCO benchmark

Pressure (bar) (mcd )UP (kg/s) (mcd )MED (kg/s) (mcd )LOW (kg/s) ave (K) T∞ Xvave (molar f.)

Experiment

CFD results

Deviation (%)

3.6 ± 0.1 36–37 32–37 30–31 128–130 0.62–0.64

3.79 35.7 30.8 33.5 135.9 0.644

5.3 −2.2 −10.7 9.8 5.3 1.6

4.1.4. MICOCO test The MICOCO benchmark (Blumenfeld et al., 2003) was conducted in a large enough facility (MISTRA) with the presence of several phenomena: plumes, natural circulation and condensation. The facility consists of a 100 m3 cylindrical vessel (4.25 m inner diameter and 7.3 m high) with thermally controlled walls and steam/helium injection lines. Three vertically stacked (upper, medium and lower) cylindrical condensers are set up inside the vessel (3.8 m diameter and 69 m2 of total surface area). The experimental procedure consists of three steps. A preheating phase, in which hot steam (195 ◦ C) is injected and the condensers are heated up to 140 ◦ C. When the atmosphere has reached 140 ◦ C and 5 bars, the steam injection is stopped and the condensers are cooled down to 120 ◦ C. In the last phase, the steam is injected at the rate of 0.12 kg/s at a dJ = 0.2 m diameter nozzle in order to reach the steady-state in about 3 h. In this paper, simulations of the last phase are shown, in order to illustrate the validation of the condensation models. A 2D cylindrical mesh (73 × 76 nodes) was used as the optimal one. The results of the main global variables have been summarised in Table 3 in comparison to the experimental data at the steady-state phase. The main variable for checking the condensation model, namely the pressure, provided good results, in the order of the experimental uncertainties. Moreover, the distribution of the condensation rate at the upper and medium condensers stayed also within the uncertainties. However, the condensation rate was over-estimated at the lower condenser being, in contrast with the experiment trends, in the order of the upper condenser and larger than the medium condenser. The reason was that the correlation is based only on the air bulk concentration and the well-mixed atmosphere conditions calculated in the experiment

Fig. 5. Terasaka’s condensation model affected by a pressure coefficient in comparison with the experimental data of the Dehbi’s tests.

The experiments performed by Anderson et al. (1998) have been used as a second database for the condensation model benchmarking. The facility is a box 2.7 m high, 1.7 m long and 0.32 m wide, with the condensers located at one upper corner. The condensers are 0.91 m long, one horizontal and the other one vertical. The experimental procedure consisted of achieving different bulk temperatures and pressures by increasing the amount of hot steam injected. A 3500-node 2D Cartesian mesh was used as the optimal one. Eleven tests were simulated (Table 2), five of them with helium in the non-condensable mixture. The Terasaka’s condensation model showed good results at the tests performed without helium. Differences with the experimental data (Table 2) lay in the order of the experimental uncertainties (±20%). No important effects from non-condensable gas concentration and pressure were identified. However, these differences increased up to 60% at higher helium mass concentrations (over to 15%). The experiment showed an accumulation of helium at the upper vessel, restricting the access of steam towards the condenser, thereby reducing the heat transfer (Anderson et al., 1998). This flow behaviour was reproduced by the CFD code, when the simulation transient was long enough in time. However, the exact conditions of the test were not well reproduced, being not homogeneous conditions: this could be the reason for the observed deviations.

Table 2 Experimental conditions, test results and comparison with the condensation model results at the Anderson’s test Test

T∞ (◦ C)

Tw (◦ C)

p (bar)

YAIR

YHe

hEXP (W/(m2 K)) t

hCFD (W/(m2 K)) t

A01 A02 A03 A05 H01 H02 H03 H04 H05 H07 H08

80 70 73.7 106.2 117.7 114.6 104.6 105.8 91.9 90 90

30 30 33.7 66.2 77.1 71.1 59.9 61.3 45.9 20 20

1 1 1.5 2.5 3 3 3 3 3 1 1

0.64 0.78 0.83 0.61 0.50 0.53 0.66 0.62 0.71 0.42 0.24

0 0 0 0 0 0.007 0.016 0.025 0.056 0 0.039

186.7 99.9 131.8 420.3 583.9 498.3 319.4 319.1 371.2 450.0 271.4

178.9 115.7 103.0 319.8 499.4 424.2 281.2 279.8 156.3 502.4 121.0

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

Fig. 6. Comparison between calculated temperature and experimental results in a vertical line located at a medium distance between the vessel axis and the condensers during the steady-state phase, at the MICOCO test.

justified this behaviour. The averaged bulk temperature was over-estimated by 6 ◦ C, a difference justified by Blumenfeld et al. (2003) on the basis of simulation of the containment walls (apart from condensers) as adiabatic walls. Although these walls were isolated with rock wool, a 15% of spurious condensation was estimated once the experiment was finished. This outcome is reflected in Fig. 6, where the experimental evolution was calculated well but over-estimated in absolute values. By contrast, the steam molar fraction at the bulk was calculated within the uncertainty range of the experimental data. This is shown in Table 3 for global values and in Fig. 7 for the local values. 4.2. NUPEC mixing tests The next step consisted on evaluating the code and the implemented models in conditions similar to the application ones. The NUPEC series of hydrogen mixing and distribution tests was selected. It is a large facility 1300 m3 in volume (14.4 m high, 10.8 m in diameter), representing a typical four-loop PWR containment at 1/4 linear scale, 25 compartments, representing those of a real containment (dome, annular rooms, loop rooms, pressuriser room, etc.) are arranged in its inside. Nozzles at dif-

639

Fig. 7. Comparison between calculated steam molar fraction and experimental results in a vertical line located at a medium distance between the vessel axis and the condensers during the steady-state phase, at the MICOCO test.

ferent locations are available to release helium (as hydrogen simulating) and steam. Moreover, 21 nozzles are located at the dome ceiling to simulate the effect of spray system in some tests. A measurement system records the transient data of atmosphere and wall temperature, helium concentration and pressure (CSNI, 1994). Nine experimental series were performed in this facility, the aim was evaluating the mixing and distribution process of gases induced by the differences on molecular weight, temperature, and other effects as natural convection, release point location and spray systems. Even though many of these tests are proprietary, five of them have been published in some way (CSNI, 1994; Stamps, 1995). The two tests without sprays were selected. They consisted of steam and helium releases at the lower and medium parts of the containment, M-4-3 and M-8-1 tests, respectively (Table 4). The very detailed nodalisation of the NUPEC containment was built up. It included all the rooms of the actual facility (except for the cavity) with a volume deviation below 0.3%. Connection pathways between compartments have been also included, trying to follow the actual locations and passage areas. For the geometrical restrictions in a structured mesh scheme, some of these connections were built up larger than the real ones

Table 4 Initial and boundary conditions for the selected NUPEC mixing tests Test

M-4-3

M-8-1

Initial pressure (kPa) Atmosphere initial temperature (◦ C) Wall initial temperature (◦ C) Injected gas Injection compartment Injected helium temperature (◦ C) Injected steam temperature (◦ C) Helium and steam injection duration (min) Helium injection rate (kg/s) Steam injection rate (kg/s) Spray operation Preheating phase

101 27–32 27–32 Helium and steam SG basement of loop D (#8) 20 110 30 0.027 0.33 No No

101 7–10 7–10 Helium and steam Pressuriser room (#22) 10 108 30 0.027 0.33 No No

640

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

Fig. 8. CFX-4 nodalisation of the NUPEC experimental facility, (left) showing the rooms, operation floors and the release points (dark squares) and (right) mesh at a vertical plane.

and the open flow was reduced further by porous media, e.g. the annular gaps. The release points were simulated in their corresponding rooms as mass flow boundary conditions. Finally, the internal and external walls were included in the geometry considering their current materials where the 1D heat conduction equation was solved. Fig. 8 outlines this geometry. A computational mesh of 220,000 nodes (characteristic cell length ranging from 0.08 to 0.5 m) in this study was refined near the release points and the connections. The first results to be compared were the pressure calculated by the code, which is an indication of the accuracy of the condensation model. For the M-4-3 test, the code reproduced very well the pressure evolution, being 3% the maximum deviation at the 400 s of the transient and being very close at the end of the experiment. For the M-8-1 test, results worsened, being the maximum deviation 6% at the end of the transient. This major deviation, respect to the previous test, is explained by the stratification conditions observed in the M-8-1 test. As observed at the validation with the Anderson’s test, the condensation model shows some deficiencies in these conditions. It is important to notice that simulations performed with MELCOR code of these experiments (Mart´ın-Valdepe˜nas et al., 2003), presented larger pressure deviations at the end of the transient, 10 and 13%, respectively. The atmosphere temperature and composition predictions showed similar behaviour results for each test. The well-mixed conditions of the M-4-3 test were suitably predicted by the CFD code, as is shown in Fig. 9. However, test measurements of the helium concentration were 15% higher than the calculated ones at the dome and lower annulus. At the release room, located at the lower part of the containment, the code reproduced the same behaviour as during the experiment.

The stratified conditions in the M-8-1 test were also predicted (Fig. 10), being the helium molar fraction deviation in the order of 11% at the dome, but much closer at the lower annulus. At the release room, located at the medium part of containment, the code did not reproduce well the experimental behaviour. The helium accumulation therein was simulated too early and, after an accumulation peak that followed the data, the concentration was established at 50% below these. Apart from the difficulties in building an exact geometry and connections for this badly connected room, some possible deficiencies in the experimental results could explain this deviation, being the atmosphere temperature at this room well estimated (Fig. 12).

Fig. 9. Comparison between calculated helium molar fraction and experimental results at the release room (#15, loop D basement), the dome (#25) and the lower annular compartment (#4), at the NUPEC M-4-3 test.

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

Fig. 10. Comparison between calculated helium molar fraction and experimental results at the release room (#22, pressuriser), the dome (#25) and the lower annular compartment (#4), at the NUPEC M-8-1 test.

641

Fig. 12. Comparison between calculated atmosphere temperature and experimental results at the release room (#22, pressuriser), the dome (#25) and the lower annular compartment (#4), at the NUPEC M-8-1 test.

Predictions for temperature behaved similar to the helium distribution: homogeneous (M-4-3) and stratified (M-8-1), see Figs. 11 and 12. Results at the dome were very close to the experimental ones in both tests, suggesting the possibility of local effects by the location of the helium sensors whose results were not so good there. By contrast, deviations were higher at the lower annulus and, especially for the M-4-3, a thermal stratification appeared in this experiment. The reason was the presence of a water pool which cannot be simulated by the CFD code, but was evidenced by complementary MELCOR calculations (Barroso, 2003). Concerning the release rooms, the code predicted temperature values much closer to the experimental ones in both tests. Finally, the performance of the 1D heat conduction model implemented in the code is shown in Fig. 13. The code reproduced well the wall temperatures of both experiments, the major Fig. 13. Comparison between calculated wall temperature and experimental results at the lower annular compartment (#6) and dome ceiling (#25), at the NUPEC M-4-3 and M-8-1 tests.

heating of the dome ceiling (hot plume impingement) and the minor of the lower containment, being these stratification more important at the M-8-1 tests by the its inhomogeneous conditions. 5. Plant application The last step is the application of the modified CFD code to actual plants. The aim of this chapter is just demonstrating and

Table 5 Containment data of both PWR plant designs

Fig. 11. Comparison between calculated atmosphere temperature and experimental results at the release room (#15, loop D basement), the dome (#25) and the lower annular compartment (#4), at the NUPEC M-4-3 test.

Free volume (m3 ) Volume over SGs (%) No. of compartments Connection SGs-dome Mass of zirconium (kg)

Westinghouse

KWU

58,000 67 ∼15 Open area 800

64,000 31 ∼60 Flaps 1020

642

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

illustrating the capacity of the code in carrying out hydrogen safety analysis in plants, rather than performing an exhaustive study in the selected plants. These studies can be reviewed in Mart´ın-Valdepe˜nas and Jim´enez (2002) and Jim´enez et al. (2003b). Two PWR plants have been selected: three-loop Westinghouse and KWU design plants. Both of them represent two different concepts of large-dry containments (Table 5): the first one is an open, lowly compartmentalised containment whereas the second one is a less open geometry with a higher degree of compartmentalisation. Moreover, the KWU one presents accidental sequences with higher hydrogen releases, being a larger mass of zirconium in the core. The sequence selected to illustrate the paper is the in-vessel phase of a SB-LOCA. The release point was located at the lower part of one loop room, particularly at the loop pump body. However, the hydrogen/steam release profiles used in this simulation were different in both plants. They were selected from results of different analysis in the reactor vessel and primary circuit compiled for each plant (Jim´enez et al., 2001, 2003b). Data for these sequences are compiled in Table 6. The main difference is the hydrogen release profile assumed for mass flow rate in each sequence: a triangular shape followed by a constant value (Westinghouse), two triangular profiles with higher and lower intense hydrogen release (KWU). Notice that the sequence simulated started with high steam concentration at the containment, meaning that the water/steam release phase of the accident was not simulated, since conditions were taken from previous MELCOR simulations, see the preceding references.

Table 6 Characteristics of the hydrogen/steam release profiles for both plants

H2 profiles Max. H2 rates (kg/s) Total H2 release (kg) Steam profiles Max. steam (kg/s) Release duration (s) Initial XH2 O (vol.%)

Westinghouse

KWU

Triangular/constant 0.913/0.023 204 Constant 4 300/2900 31

Triangular/triangular 0.60/0.71 797 Constant/triangular 5.32/5.54 2000/500 45

5.1. Geometry and results at the Westinghouse design plant A 3D geometry was built with the main rooms, free volumes and flow junctions included. However, the metal-made structures (SGs, pressuriser, loop ducts, etc.) were removed for saving cell size and geometry resources. The characteristic cell size was approximately 1 m3 , amounting a total of 90,000 cells, Fig. 14. The nine main compartments (the three loops, dome, annulus, pressuriser, pressuriser relief tank, refuelling reactor pool, sump) were selected via USER3D facility to apply the combustion criteria. Results in Fig. 15 show the way hydrogen and steam flow – initially released as a jet through the break towards the dome – loosing their inertia soon within the break room and travelled upward as a plume driven by buoyancy. A poorly mixed distribution map of hydrogen and steam was established at this room, with the richest area at the upper part and very few hydrogen underneath the release point. Once the plume entered the dome-free region, it did not mix within the full area, being

Fig. 14. Geometry and grid for the CFD analysis of the PWR Westinghouse design plant case: (left) perspective of the containment free volume, (right) inner compartments: the three loops, pressuriser and relief tank rooms and refuelling reactor pool.

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

643

Fig. 15. Fields of hydrogen molar fraction at a vertical plane through the break room of the PWR Westinghouse case. Temporal evolution showing the plume formation, the plume spreading at the dome ceiling, the mixing at the dome and the final homogeneous condition. (Notice that saturated black correspond with the concentrations around and over 10 vol.%, jet close to the release point.)

the plume buoyancy forces larger than the natural convection patterns. Subsequently, the plume impinged the containment ceiling and spread out. After the hydrogen release maximum, due to the low elevation of the release point and the continuous release of steam, the hydrogen was driven to the lower part of the containment, leading to homogeneous conditions. Application of the criteria showed flammable mixtures only at the release room, at the dome during the release peak and the whole containment at the end of the constant release phase. The maximum accumulations of hydrogen happened at the break room and dome during the release peak. Neither FA nor DDT was possible in any compartment during the analysis, since the GLOBAL application of criteria showed values much below the threshold, being iσ < 0.75 and iλ < 0.1. However, application of the LOCAL criteria at the release room showed some clouds

which fulfilled the criteria during some few minutes around the time of maximum release and involving clouds of some 20 kg of hydrogen (Fig. 16). 5.2. Geometry and results at the KWU design plant The 3D geometry included the main rooms, free volumes, connections and the main metal-made structures (SGs, main pumps, pressuriser and relief tank, hot and cold legs and accumulators), see Fig. 17. A computational mesh of 214,000 nodes (0.003–9 m3 cell sizes) was applied in this study, with finest meshing at the release rooms, near the metal-made equipment and the release point. Twelve main compartments were selected as USER3D regions to apply the combustion criteria (the three loops, dome, annulus, pressuriser, pressuriser

644

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

Fig. 16. Temporal evolution of the FA and DDT index at the break room during the peak phase and first seconds of the constant phase in the PWR Westinghouse case.

relief valves, refuelling reactor pool, sump, the three staircases and auxiliary room). Notice that, in terms of simplification of the criterion application, the staircase and auxiliary regions lump several rooms together. The behaviour of the hydrogen in this enclosure was similar to the Westinghouse case: the light gas rose up as buoyant plume, spreading at the dome and finally reaching homogeneous conditions, but gathering in the short term at the release room. The main differences show at the second release peak, when the hydrogen plume is mixed in an atmosphere with 7 vol.% of hydrogen, Fig. 18 illustrates this behaviour.

Results of this case are summarised in Table 7. During the first peak at the release region (break loop), the combustion indices stayed below the threshold values, except for a small cloud of approximately 1 kg on DDT LOCAL conditions. The rest of the rooms showed values similar to that of the dome, i.e. always below the threshold for both indices. During the second peak, with higher mixture sensitivity, the GLOBAL and LOCAL values of iσ reached up to 0.9 at the dome and were in fact above the threshold at the release room. The GLOBAL DDT was discarded at any room but a small cloud reached DDT conditions near the break point.

Fig. 17. Geometry for the CFD analysis of the PWR KWU design plant: (left) metallic equipments and release point and (right) concrete structures seeing the operation floor and biological shield.

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

645

Fig. 18. Fields of hydrogen molar fraction at a vertical plane through the break room of the PWR KWU case. Temporal evolution showing the plume spreading at the dome ceiling and the final homogeneous condition at both release peaks.

Table 7 LOCAL and GLOBAL results of the FA and DDT index and the more relevant regions of the PWR KWU containment USER3D

Maximum FA index

Maximum DDT index

kg H2

GLOBAL

LOCAL

GLOBAL

LOCAL

First hydrogen peak Break loop Dome

0.76 0.76

0.85 0.76

0.002 0

2.97 0

22.2/0.8 18.2/0

Second Hydrogen peak Break loop Dome

0.92 0.90

0.92 0.90

0.027 0.035

16.3 0.12

35.6/1.5 273/0.03

kg H2 , hydrogen mass within the region at flammable conditions/at the cloud with maximum DDT index.

6. Conclusions A CFD code for analysis of hydrogen behaviour within NPP containments during a severe accident is presented in this paper based on the commercial multi-purpose code, CFX-4. Because of the lacks in simulating important phenomena in the containment some models have been implemented in the code. The main are: (1) steam condensation onto walls in presence of non-condensable gases, (2) heat conduction, (3) fog and rain

formation, (4) material properties and (5) hydrogen combustion criteria. Selected validation exercises were undertaken next. They have demonstrated the capacity of the code in simulating the most relevant phenomena in the hydrogen behaviour. The basic experiments simulated have shown the code abilities in predicting the stratification phenomena and the mixing process at relatively small enclosures. The condensation phenomena were predicted properly when the Terasaka and Makita model

646

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647

was applied. Additionally, the condensation benchmarks have permitted to identify the worst working conditions for the model: high helium (hydrogen simulating) concentrations and high absolute pressures. The second effect has been overcome by a correction factor p2/3 , which is in the order of the condensation dependence with pressure. The code benchmarking against the MICOCO test verified the capacity of the code in simulating steam plumes and condensation in larger facilities. A second step in the validation procedure consisted of the simulation of the integral mixing and distribution M-4-3 and M8-1 tests within the NUPEC large-scale facility, without sprays and the code reproduced suitably the opposite conditions of both tests: homogeneous in the M-4-3 and stratified in the M-8-1. Moreover, the code predictions for pressure, helium concentration, atmosphere temperature and wall temperature evolutions lay very close to the experimental results. Only at the badly connected release room during the M-8-1 test, the helium concentration was extremely underestimated, although this is not a fully convincing experimental result. Two plant applications have exemplified the use of the code in two kinds of PWR containment designs: Westinghouse and KWU, as representative of two different concepts in large-dry containments. The first one is an open-junction, lowly compartmentalised containment whereas the second one is a less open geometry with a higher degree of compartmentalisation. Even though the sequences simulated were not fully comparable to each other, this analysis shows an illustrative exercise. Both designs featured local hydrogen clouds at the release rooms and possible FA or DDT conditions during the maximum hydrogen release rates. Moreover, at the KWU design worse conditions took place at the whole containment, especially if a second intensive zirconium oxidation (e.g. by core reflooding) would happen. Finally, new tasks would be recommended as future works. They would involve the three main aspects in this paper: model implementation, validation and plant application by models simulating water spray performance, pool formation and PAR actuation, and the extension of the database for carbon monoxide and carbon dioxide. Moreover, further model developments in condensation and fog or rain formation would also be of concern, being validated against representative tests, especially integral ones. Systematic code application to both kinds of plant during representative accident sequences would allow evaluating the containment behaviour at the expected conditions and comparing both design performances against the hydrogen problem. Acknowledgments Authors wish to acknowledge to the Chair of Nuclear Technology, Universidad Polit´ecnica de Madrid, where this research work was performed. Moreover, acknowledgement is extended to the support of the Fifth Framework Program of the EU under the EXPRO project for the NUPEC validation and KWU plant application, and the Spanish Nuclear Safety Council (CSN) for the Westinghouse plant application.

References AEA, 2001. CFX 4.4: Solver. CD-ROM, CFX International, AEA Technology plc., Harwell, United Kingdom. Anderson, M.H., Herranz, L.E., Corradini, M.L., 1998. Experimental analysis of heat transfer within AP600 containment under postulated accident conditions. Nucl. Eng. Des. 185, 153–172. Asano, K., Nakano, Y., Inaba, M., 1979. Forced convection film condensation of vapors in the presence of non-condensable gas on a small vertical flat plate. J. Chem. Eng. Jpn 12 (3), 196–202. Bird, B.R., Stewart, W.E., Lightfoot, E.N., 1960. Transport Phenomena. John Wiley & Sons, New York. Blumenfeld, L., Paill`ere, H., Choi, Y.J., Studer, E., Baumann, W., Travis, J.R., Chin, Y.S., Krause, M., George, T.L., Wheeler, C.L., Mart´ın-Valdepe˜nas, J.M., Putz, F., Schwarz, S., 2003. CFD-simulation of mixed convection and condensation in a reactor containment/the MICOCO benchmark. In: Proceedings of the 10th International Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-10), Seoul. Breitung, W., 1997. The analysis of hydrogen behaviour in severe accidents. In: Alonso, A. (Ed.), EUROCOURSE-97: Analysis of Severe Accidents in Light Water Reactors, Lecture 11, F2 I2 , EC-Nuclear Fission Safety Programme (1994–1998). Breitung, W., Royl, P., 2000. Procedure and tools for deterministic analysis and control of hydrogen behaviour in severe accidents. Nucl. Eng. Des. 202, 249–268. Broughton, J.M., Kuan, P., Petti, D.A., Tolman, E.L., 1989. A scenario of the Three Mile Island Unit 2 accident. Nucl. Technol. 87, 34–53. Byun, C.S., Jerng, D.W., Todreas, N.E., Driscoll, M.J., 2000. Conceptual design and analysis of a semi-passive containment cooling system for a large concrete containment. Nucl. Eng. Des. 199, 227–242. CAB, 1988. Programa TABAGUA: Tabla de Propiedades del Agua Liviana. Centro At´omico de Bariloche. Chapman, A.J., 1984. Transmisi´on del Calor, third ed. Bellisco, Madrid, pp. 410–426. Collier, J.G., 1972. Convective Boling and Condensation. McGraw-Hill, UK, pp. 314–359. CSNI, 1994. Final Comparison Report on ISP-35: NUPEC H2 Mixing and Distribution Test (Test M-7-1). CSNI Report, OECD Nuclear Energy Agency, Nuclear Safety, NEA/CSNI/R(94)29, OECD/GD(95)29. CSNI, 1999. State of Art Report on Containment Thermalhydraulics and Hydrogen Distribution. CSNI Report, OECD Nuclear Energy Agency, Nuclear Safety, NEA/CSNI/R(99)16. CSNI, 2000. State of Art Report on Flame Acceleration and Deflagration to Detonation Transition in Nuclear Safety. CSNI Report, OECD Nuclear Energy Agency, Nuclear Safety, NEA/CSNI/R(2000)7. Dehbi, A.A., Golay, M.W., Kazimi, M.S., 1991. Condensation experiments in steam–air, steam–air helium mixtures under turbulent natural convection. In: National Conference of Heat Transfer, AIChE Symposium, pp. 19–28. Dorofeev, S.B., Kotchourko, A.S., Chaivanob, B.B., 1989. Detonation onset conditions in spatially non-uniform combustible mixtures. In: Proceedings of Sixth International Symposium on Loss Prevention and Safety Promotion in the Process Industries, vol. 4, Oslo, pp. 221– 229. Dorofeev, S.B., Efimenko, A.A., Bezmelnitsin, A.V., 1996. Analysis and Evaluation of DDT Criteria for Severe Accidents Conditions. KIM Report, Final report for FZK. Dorofeev, S.B., Kuznetsov, M.S., Alekseev, V.I., Efimenko, A.A., Breitung, W., 2001. Evaluation of limits for effective flame acceleration in hydrogen mixtures. J. Loss Prevent. Process Ind. 14, 583–589. EUR, 2003. Hydrogen Energy and Fuel Cells: A Vision of Our Future. Final Report of High Level Group, EUR 20719 EN. Gallego, E., et al., 2005. An intercomparison exercise on the capabilities of CFD models to predict distribution and mixing of H2 in a closed vessel. First International Conference on Hydrogen Safety, HySafe Network of Excellence, Pisa. Gauntt, R., et al., 2000. MELCOR Computer Code Manuals. Sandia NL Report, NUREG/CR-6119.

J.M. Mart´ın-Valdepe˜nas et al. / Nuclear Engineering and Design 237 (2007) 627–647 Gordon, S., McBride, B.J., 1994. Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications I. Analysis. NASA Report, RP-1311. Henrie, J.O., Postma, A.K., 1983. Analysis of the Three Mile Island Unit Hydrogen Burn. Idaho NL Report, GEND-INF-023-vol. 4, DE83 012186. Henrie, J.O., Postma, A.K., 1987. Lessons Learned from Hydrogen Generation and Burning During the TMI-2 event. Idaho NL Report, GEND-061, DE87 010696. Herranz, L.E., Anderson, M.H., Corradini, M.L., 1998. Diffusion layer model for steam condensation within the AP600 containment. Nucl. Eng. Des. 183, 133–150. Incropera, F.P., DeWitt, D., 1996. Fundamentals of Heat and Mass Transfer, fourth ed. John Wiley & Sons, New York, pp. 556–564. ITER, 2001. Generic Site Safety Report (GSSR): Ultimate Safety Margin. ITER Report, GSSR, G 84 RI 7 01-07-10 R 1.0. Jim´enez, M.A., Mart´ın-Valdepe˜nas, J.M., Fern´andez, J.A., Mart´ın-Fuertes, F., 2001. An´alisis con MELCOR 1.8.4 del Riesgo de Combusti´on de H2 en la Contenci´on de CN Asc´o para la Envolvente de Secuencias SBLOCA. CTN/UPM Report, CTN-13/01 (in Spanish). Jim´enez, M.A., Mart´ın-Valdepe˜nas, J.M., Coe, I., Pears, K., 2003a. A Methodology for Hydrogen Combustion Risk Assessment. Deliverable D27 EXPRO Project (V EU-FWP), CONTRACT EVG1-CT-2001-00042. Jim´enez, M.A., Mart´ın-Valdepe˜nas, J.M., Coe, I., Peers, K., 2003b. Combustion Risk Application to EU Power Plant Enclosures. Part I Deliverable D28 EXPRO Project (V EU-FWP), CONTRACT EVG1-CT-2001-00042. Kim, M.H., Corradini, M.L., 1990. Modelling of condensation heat transfer in a reactor containment. Nucl. Eng. Des. 118, 193–212. Kumar, R.K., 1985. Flammability limits of hydrogen–oxygen–diluent mixtures. J. Fire Sci. 3 (4), 245–262. Marshal, B.W., 1986. Hydrogen:Air:Steam Flammability Limits and Combustion Characteristics in the FITS Vessel. Sandia NL Report, NUREG/CR3468, SAND84-0383. Mart´ın-Valdepe˜nas, J.M., Jim´enez, M.A., 2002. CFD analyses of hydrogen risk within PWR containments. In: Proceedings of International Technical Meeting on Use of CFD Codes for Safety Analyses of Reactor Systems, Including Containment, IAEA-OECD/NEA, Pisa. Mart´ın-Valdepe˜nas, J.M., Jim´enez, M.A., Barroso, R., Coe, I.M., Peers, K.L., Rowney, B.A., 2003. MELCOR, CFX-4, COMPACT Benchmarking Exercise on NUPEC Mixing Test. Part II Deliverable D28 EXPRO Project (V EU-FWP), CONTRACT EVG1-CT-2001-00042. Mart´ın-Valdepe˜nas, J.M., 2004. Numerical Models Implemented in a Computational Fluid-Dynamics Code for Hydorgen Combusti´on Risk Assessment within Containments. Ph.D. Thesis, Universidad Polit´ecnica de Madrid, Madrid, ISBN 84-609-3478-0 (in Spanish). Mart´ın-Valdepe˜nas, J.M., Jim´enez, M.A., Mart´ın-Fuertes, F., Fern´andez, J.A., 2005. Comparison of film condensation models in presence of noncondensable gases implemented in a CFD code. Heat Mass Transfer 41 (11), 961–976. Paill`ere, H., et al., 2003. Simulation of H2 release and combustion in large scale geometries: models and methods. In: Proceedings of Super-Computing for Nuclear Applications, SNA-2003, Paris.

647

Peterson, P.F., 1996. Theoretical basis for the Uchida correlation for condensation in reactor containments. Nucl. Eng. Des. 162, 301–306. PHEBUS, 1997. FPT0 Final Report. IPSN Report, Cadarache, IPSN/DRS/SEA/ LERES/97/1815 CLT.47 12 40. Reid, R.C., Prausnitz, J.M., Sherwood, T.K., 1988. The Properties of Gases and Liquids. McGraw Hill, New York. Renson, C., Goldstein, S., 1987. Stratification d’Helium dans une Enceinte. Premiers Resultats des Essais Isothermes, CEA Report, S´aclay, Rapport DEMT/87/127. Royl, P., Travis, J.R., Haytcher, E.A., Wilkening, H., 1996. Analysis of mitigation measures during steam/hydrogen distributions in nuclear reactor containments with the 3D field code GASFLOW. In: Proceedings of OECD/NEA/CSNI Workshop on the Implementation of Hydrogen Mitigation Techniques, AECL-11762, NEA/CSNI/R(96)8, Winnipeg. Royl, P., Travis, J.R., Breitung, W., Seyffarth, L., 1997. Simulation of hydrogen transport with mitigation using the 3D field code GASFLOW. In: Proceedings of International Meeting on Advanced Reactors Safety, Orlando, ISBN 0-89448-624-1. Royl, P., Rochholz, H., Travis, J.R., Necker, G., 1999. Three-dimensional analysis of steam/hydrogen transport with catalytic recombiners in nuclear reactor containments using the computer code GASFLOW. In: Proceedings CSARP Meeting, Albuquerque. Royl, P., Rochholz, H., Breitung, W., Travis, J.R., Necker, G., 2000. Analysis of steam and hydrogen distributions with PAR mitigation in NPP containments. Nucl. Eng. Des. 202, 231–248. Royl, P., Travis, J.R., Breitung, W., 2001. Modelling of catalytic and application in 3D containment using the 3D code GASFLOW. In: Proceedings of IAEA Technical Committee Meeting on Implementation of Hydrogen Techniques and Filtered Containment Venting, Cologne. Royl, P., Travis, J.R., Baumann, W., Breitung, W., et al., 2002. Status of development, validation, and application of the 3D CFD code GASFLOW at FZK. In: Proceedings of IAEA-OECDE Technical Meeting. Shebeko, Y.N., Keller, V.D., Yeremenko, O.Y., Smolin, I.M., Serkin, M.A., Korolchenko, A.Y., 1988. Regularities of formation and combustion of local hydrogen–air mixtures in a large volume. Chem. Ind. 21, 24–27 (in Russian). Stamps, D.W., Berman, M., 1991. High-temperature hydrogen combustion in reactor safety applications. Nucl. Sci. Eng. 109, 39– 48. Stamps, D.W., 1995. CONTAIN Assessment of the NUPEC Mixing Experiments. Sandia NL Report, SAND94-2880 UC-505. Terasaka, H., Makita, A., 1997. Numerical analysis of the PHEBUS containment thermal hydraulics. J. Nucl. Sci. Technol. 34 (7), 666– 678. Travis, J.R., Spore, J.W., Royl, P., et al., 1998. GASFLOW: A Computational Fluid Dynamics Code for Gases, Aerosols, and Combustion. Los Alamos and FZK Report, vols. 1–3, LA-13357-M, FZKA-5994. Vierow, K.M., Schrock, V.E., 1991. Condensation in a natural circulation loop with noncondensable gases. Part-1: heat transfer. In: Proceedings of the International Conference on Multiphase Flows’91, Tasukaba, pp. 183– 186.