Scholarly article on topic 'Validation of NEPTUNE-CFD Two-Phase Flow Models Using Experimental Data'

Validation of NEPTUNE-CFD Two-Phase Flow Models Using Experimental Data Academic research paper on "Chemical engineering"

CC BY
0
0
Share paper
OECD Field of science
Keywords
{""}

Academic research paper on topic "Validation of NEPTUNE-CFD Two-Phase Flow Models Using Experimental Data"

Hindawi Publishing Corporation

Science and Technology of Nuclear Installations

Volume 2014, Article ID 185950, 19 pages

http://dx.doi.org/10.1155/2014/185950

Research Article

Validation of NEPTUNE-CFD Two-Phase Flow Models Using Experimental Data

Jorge Pérez Mañes,1 Victor Hugo Sánchez Espinoza,2 Sergio Chiva Vicent,3 Michael Böttcher,2 and Robert Stieglitz2

1 Laboratoire d'Etudes et de Simulation des Systémes, CEA Cadarache, CAD/DEN/DER/SESI, Bat. 212, 13108 St. Paul Lez Durance Cedex, France

2 Karlsruhe Institute of Technology, Institute for Neutron Physics and Reactor Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

3 Mechanical Engineering and Construction Department, Jaume I University, Avenida de Vicent Sos Baynat, s/n, 12071 Castellon, Spain

Correspondence should be addressed to Victor Hugo Sánchez Espinoza; victor.sanchez@kit.edu Received 16 February 2014; Accepted 9 April 2014; Published 30 June 2014 Academic Editor: Eugenijus Uspuras

Copyright © 2014 Jorge Perez Manes et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This paper deals with the validation of the two-phase flow models of the CFD code NEPTUNEC-CFD using experimental data provided by the OECD BWR BFBT and PSBT Benchmark. Since the two-phase models of CFD codes are extensively being improved, the validation is a key step for the acceptability of such codes. The validation work is performed in the frame of the European NURISP Project and it was focused on the steady state and transient void fraction tests. The influence of different NEPTUNE-CFD model parameters on the void fraction prediction is investigated and discussed in detail. Due to the coupling of heat conduction solver SYRTHES with NEPTUNE-CFD, the description of the coupled fluid dynamics and heat transfer between the fuel rod and the fluid is improved significantly. The averaged void fraction predicted by NEPTUNE-CFD for selected PSBT and BFBT tests is in good agreement with the experimental data. Finally, areas for future improvements of the NEPTUNE-CFD code were identified, too.

1. Introduction

The validation of the two-phase flow modelling capability of CFD codes, for example, NEPTUNE-CFD, is mandatory for its application in the design and safety evaluation of energy systems. The goal thereby is to demonstrate that the two-phase flow models of CFD codes are able to predict the most relevant flow regimes under pre- and postcritical heat flux (CHF) conditions. The focus is on the accurate prediction of the pressure drop, void fraction, critical power, departure from nucleated boiling (DNB), and so forth. The ways how the CFD codes are modelling the heat transfer between a solid and the coolant (wall heat transfer) and the liquid-vapour interphase heat transfer (bulk heat transfer) differ from code to code. The validation of the two-phase flow heat transfer models of NEPTUNE-CFD requires the

understanding of the implemented mathematical-physical models in the code as well as their interaction with the heat conduction models of the solids and the fluid dynamics. In the medium term a combined application of CFD and system or subchannel codes will lead to a more realistic prediction of safety-relevant phenomena in nuclear reactors. The simulations performed during this work are focused on the void fraction prediction in rod bundles of light water reactors (LWR) using the NEPTUNE-CFD 1.0.8. In this paper, the investigations performed to validate the two-phase flow models of NEPTUNE-CFD under steady-state and transient conditions using the database provided by the PSBT and BFBT tests are presented and discussed. This research code is being developed and tested within the European NURISP and NURESAFE projects. First of all, the main features of NEPTUNE-CFD are presented in detail. Then, the

validation work based on both PSBT and BFBT is discussed in detail. Finally, a summary and the main conclusions are given.

2. The NEPTUNE-CFD Code

The NEPTUNE-CFD solver is based on a pressure correction approach to simulate multicomponent multiphase flows by solving a set of three balance equations for each field (fluid and/or gas phase) in a Reynolds Averaged Navier Stokes (RANS) approach. These fields can represent many kinds of multiphase flows; among them also is bubbly flow. The solver is based on a finite volume discretization together with a collocated arrangement for all variables. The two-fluid models of the NEPTUNE-CFD are designed specifically for the simulation of two-phase transients in nuclear reactors. The models and closure laws used by NEPTUNE-CFD [1] presented in this paper are focused only on the description of the boiling phenomena and heat flux partitioning at the wall under boiling conditions.

When an averaging operation is performed, the major part of the local information at the interfaces and the physics governing the different types of exchanges at a microscale are lost. As a consequence a number of closure relations (also called constitutive relations) must be supplied to close the balance equations so that they can be mathematically solved. Three different types of closure relations can be distinguished: interfacial mass and heat transfer terms (i.e., the molecular and turbulent transfer terms) and the wall heat transfer terms. In the next paragraphs, a short description of the interfacial and wall transfer terms important for the description for the boiling phenomena will be presented. Details about the NEPTUNE-CFD models can be found in the theory and user manual [2, 3].

2.1. Heat and Mass Transfer. If the mechanical terms are neglected in comparison to the thermal terms in the averaged form of the energy jump condition, this condition reduces to

• Hti +

<?"; • A¡)

interfacial heat transfer term is the product of the interfacial area concentration by the interfacial heat flux density ;):

<?"; = hk¡ • (Xat (P) - Tk) > K =

The interfacial heat flux density can be defined as (3), where hki denotes a heat transfer coefficient, Tk the average temperature of phase k, and Tsat(p) the saturation temperature. The interfacial area concentration is expressed as A; = 6a/d, where a is the void fraction, d is the Sauter mean bubble diameter (SMD), and is the thermal conductivity of the liquid. Depending on the Jakob number given by (4), there are two possible scenarios: condensation (Ja < 0) or evaporation (Ja > 0):

Pi • Cpl • (Tl - Tsat)

Pa • L .

The thermal capacity of the liquid is Cpi, and L is the latent heat of vaporization. In case of condensation, the Nusselt number is

Nu = 2 + 0.6 • Re1/2Pr1/3, Re =

El^I Pr = t±

where Refc is the bubble Reynolds number, Pr is the liquid Prandtl number, is the liquid kinematic viscosity, and is the thermal liquid diffusivity. In case of evaporation the Nusselt number is defined as follows:

Nu = max (Nu1, Nu2, Nu3),

4 • Pe

Nu3 = 2,

where Pe is the Péclet number and it is defined as follows:

Pe = Re •Pr =

d,\V„ - V,

The heat transfer term between the vapour and the interface for the case of bubbles is written as

This important relation (together with the mass jump condition: Tc = T^) allows computing the mass transfer terms as functions of the interfacial heat transfer terms q'^ - A; and the interfacial-averaged enthalpies Hki:

rc _ rC _ % + *ïvi A

r = V = H W ^

-H,,; - -Hi,

Since no information about the dependence of the interfacial-averaged enthalpies Hki is known, two basic assumptions can be made: either the interfacial-averaged enthalpies Hki are identified to the phase-averaged ones H' or the interfacial-averaged enthalpies Hki aregiven by thesaturationenthalpies. In NEPTUNE-CDF the first assumption is made. Each

' ' PV ^ frr rp N

<?v; = av-(Tsat - Tv),

where Cpv is the gas heat capacity at constant pressure and x is a characteristic time given by the users. This relation simply ensures that the vapour temperature TV remains very close to the saturation temperature Tsat, which is the expected result for bubbly flows with sufficiently small bubbles (flow in a PWR core in conditions close to nominal ones).

2.2. Interfacial Momentum Transfer. The interfacial transfer of momentum M k appears in the balance momentum equation as a source term. It is assumed tobe the sum of five forces:

M k = MDk +MlM + MLk + MTkD + .

The five terms are the averaged drag, added mass, lift force, turbulent dispersion force, and wall lubrication force per unit volume. The drag force definition used in the case of bubbly flow is the one given by the Ishii and Zuber correlation [4], where the calculation of the drag coefficient is based on the local flow regime. The added-mass force takes into account the effect of the bubbles concentration according to Zuber [5] and Ishii [6]. The lift force describes the particular case of a weakly rotational flow around a spherical bubble in the limit of infinite Reynolds number according to Auton [7]. It has been empirically modelled by Tomiyama et al. [8]. The turbulent dispersion force tends to move the bubbles to locations with less void density; it is correlated by Lance and Lopez de Bertodano [9].

2.3. Interfacial Area Modelling. The interfacial area concentration (IAC) transport equation is given by (10) according to Yao and Morel [10]. This formulation is based on the model of Wu et al. [11]:

dt v ' -V' 3

T., - a

Dv ■ PV

12 ■ n ■ )2 ■ +

nuc j2

„ ■ w d.

The advantage of (10) is that the nucleation phenomenon clearly appears as a single term given by the product of the bubble number density source term 0nuc with the surface of a nucleated bubble. If for this particular term the SMD (d) is replaced by the bubble detachment diameter dd (for wall nucleation), the IAC is changed accordingly due to the fact that the newly nucleated bubbles often have a smaller diameter dd than the SMD. The bubble detachment diameter is calculated using Unal's correlation ((17) and (15)). Closure relations must be proposed for the bubble number density source terms 0nuc, coalescence 0coa, and breakup 0brk. An example of such models is proposed by Yao and Morel [10].

2.4. Wall Boiling Model. The nucleate boiling term (q^k) appears in (11) as a wall-to-fluid heat transfer term per unit volume and time. It is assumed that all applied heat is transferred to the liquid phase. Hence the contribution to the vapour (q'W'V) is zero. To model the wall-to-fluid heat transfer at nucleate boiling, a two-step approach is used. The two steps include the calculation of the condition for boiling incipience in terms of critical wall superheat according to the Hsu criterion [12] and calculation of heat flux partitioning. Following the analysis of Kurul and Podowski [13], the wall heat flux is split into three terms:

(i) a single phase flow convective heat flux qc at the fraction of the wall area unaffected by the presence of bubbles;

(ii) a quenching heat flux q where bubbles departure bring cold water in contact with the wall periodically;

(iii) evaporation heat flux qe needed to generate the vapour phase.

Each of these three phenomena is expressed by a heat flux density (per unit surface of the heated wall) which is related to the volumetric heat flux by the following relation:

m A w n A w / 1wl =—■ Iwl = —■ (ic + % + 1e

where Aw is the heated wall surface in contact to the cell having volume V; therefore, q'w is expressed in W/m3 and q'W as well as qc, q , and qe are expressed in W/m2. The quantities qc, qq, and qe denote the heat flux densities due to liquid convective heat transfer, quenching, and evaporation. The liquid convective heat transfer per unit surface of the heated wall is written as

Vc = Ac^ Kg ■ (Tw - Ti) >

where Tw is the wall temperature and hlog is a heat exchange coefficient which is given by

fyog - Pi ■ Cpl ■ —,

in which u* is the wall friction velocity and T+ is the nondimensional liquid temperature. The velocity u* is calculated from the logarithmic law of liquid velocity in the wall boundary layer. The nondimensional temperature follows a similar logarithmic profile.

The heat flux density due to quenching is written as

-A b■ ^q S ■

2 ■ h ■ (Tw-T¡)

n ■ av tq

where Ab is the wall fraction occupied by bubble nucleation, f is the bubble detachment frequency, tq is the quenching time, and al is the liquid thermal diffusivity. The two fractions A c and Ab are given by

A h - min n

Ar -1-A.

Here n is the active nucleation sites density (per unit surface of the heated wall) and dd is the bubble detachment diameter. The active nucleation sites density is modelled according to Kurul and Podowski [13],

- (210 ■ (Tw - Tj)

as a function of the wall superheating. The bubble detachment diameter is given by the correlation from Unal [14]. This correlation is valid for subcooled liquid but has been extended to saturated liquid. The bubble detachment diameter is given by

dd - 2.4210 5 ■ p0

where p is the pressure and a, b, and f are given by the following relations:

(TW ^sat) ■ ^s 2 ■ pv ■ L ■

Here As and as denote the wall conductivity and thermal diffusivity, respectively, pv specifies the vapour density, and L is the latent heat of vaporization. In the modified correlation b is given by

(Tsat -Tt)

2- (1 -pjpi)' 1

2- (l-pv/p¡) 0.0065 - pi - Cpl-

1c + qq + <?e

St < 0.0065 St > 0.0065, (19)

where ||V^|| is the norm of the liquid velocity and St is the Stanton number which is defined by

<?c + 1q + <?e

St = -

Pi • Cpl -||VlII- (Tsat - T;)' The quantity <p appearing in (17) is given by

<p = max (1

V0 = 0.61 m/s. (21)

The quenching time and the bubble detachment frequency are modelled as

u = 7> / =

4 - ff-

3 Pr dd '

The third heat flux density qe used for evaporation is given by

<?e = / -

rc- dd3 6

To ensure a grid independent solution, the liquid temperature Tl in the wall boiling equations is calculated from the logarithmic temperature profile in a given nondimensional distance from the wall at y+ = 250. This solution is proposed by Egorov and Menter [15]. The reason is that at the centre of the wall-adjacent cell high temperature and void fraction gradients are expected which are strongly dependent on the nodalisation. Taking into account the self-similarity, the non-dimensional temperature profile in the wall boundary layer T; at y+ = 250 reads

T; (/ = 250)

T+ (/ = 250) . + ,, (24)

= -r+((;+ = wc) - -ri-= WC))'

where the subscript WC denotes the wall-adjacent cell. This approach is valid only if the wall-adjacent cells remain in log region of the wall boundary layer (30 < y+ < 300).

2.5. Boiling Model Extension for DNB Modelling. The basic wall heat flux partitioning model assumes that the amount of water on the wall is sufficient to remove heat from the wall to be used for evaporation. Superheating of the vapour that occurs at high void fractions is not modelled. Under these assumptions, the basic heat flux partitioning model cannot

be used for CHF conditions. In order to take into account the phenomenon of temperature excursion at DNB conditions, the heat flux partitioning model can be generalized as follows:

«wall = /«1 - (<?/ + + + (1- fal) - <iv (25)

The fourth part of the wall heat flux, qv, is the diffusive heat flux used to preheat the gas phase:

ív = ^wf,v - (^wall - Tv),

where fowfV is the wall heat transfer coefficient calculated from the temperature wall function for the vapour phase and Tv is the vapour temperature at the centre of the wall-adjacent cell. fa1 is the phenomenological function, which depends on the liquid volume fraction a1 and takes care for the numerically smooth transition between nucleate boiling regime and CHF regime. The generalized model assumes function fa1 in the following form:

«1 > a1,crit : /«1 = 1 - 2 - e

«1 < a1,crit : /«1 = 1 -( )

2 V a1,crit /

-2Q(«1-«1,crit)

The extension of the wall-heat-flux-partitioning model was used to take into account the CHF condition. The local void fraction equal to 0.8 can be used as a criterion for the CHF occurrence. This value is close to the Weisman and Pei DNB criterion with the void fraction equal to 0.82 according to [16].

3. Contribution to the NEPTUNE-CFD Validation

The NEPTUNE-CFD code with the two-phase heat transfer models explained in Section 2 is used for the posttest analysis of the PWR and BWR-relevant bundle tests: PSBT and BFBT. The PSBT [17] and BFBT [18] provide void fraction measurements for conditions that are representative for LWR.

3.1. OECD/NRC NUPEC PWR PSBT Benchmark

3.1.1. Scope and Description of the Benchmark. The test bench of the PSBT experiment is shown in Figure 1. The geometry described is a single centred isolated subchannel with a pitch of 8 mm between the walls of the rods. The subchannel has 4 electrically heated walls with a heated axial length of 1.55 m and a uniform axial power distribution. Averaged void fraction data is provided over a cross-sectional area located at 1.4 m distance from the bottom of the heated section. The experimental data are collected by X-ray densitometer. For the comparison, six different experiments have been selected from the PSBT database for the simulation with NEPTUNE-CFD. These experiments belong to Exercise 1 (steady-state single subchannel benchmark) from Phase 1 (void distribution Benchmark). The main parameters of these tests such as power, inlet temperature, pressure, mass flow, and subcooled liquid temperature are specified in Table 1. The NEPTUNE-CFD code lacks a steady-state algorithm; thus a null transient is performed for these simulations.

l.crit

Coolant outlet ti

Z = Z =

1.55 m

Z = 0 m

Upper electrode Subchannel model - Measure section

Red color indicates the location of the electrically heated walls

6.07 6

Figure 1: (a) Experimental setup of the PSBT benchmark; (b) cross-sectionalcut illustrating the location of the heat source; and (c) numerical model in NEPTUNE-CFD showing the velocity field for case 1.2211.

3.1.2. Applied NEPTUNE-CFD Models. The turbulent transfer term applied is the k-e model for the liquid phase and the local equilibrium turbulence model [2] for the dispersed phase. A standard wall function adapted for two-phase flow is applied. The steam/water properties are based on the International Association for the Properties of Water and Steam (IAPWS) data [19]. The steam is set at saturation temperature, and it has a slip condition at the wall. The drag and nondrag forces explained in the Section 2.2 are applied to the simulation. The wall lubrication force is not estimated. In case of boiling at a wall, the model can consider an extra artificial quenching flux. Superficial heat flux has been imposed as a boundary condition at the heated walls. As initial conditions, a water mass flow rate at the inlet and a constant pressure at the outlet are imposed.

3.1.3. Space Discretization of the Studied Domain. The cross section of the subchannel geometry is illustrated in Figure 2. The design ofthe space discretization was performed with the tool ICEM, which is included as a part of the ANSYS-CFX package. The cells in the near wall region are thinner in order to describe the velocity and temperature gradients accurately. At those locations where the velocity is not expected to have steep gradients, the grid is coarser, for example, at the centre of the subchannel.

To catch the physical phenomena near the wall, for example, the void fraction by the numerical codes, a more refined discretization is necessary. Hence, four different spatial resolutions of the subchannel were tested, all of them consist of structured meshes in order to avoid diffusivity problems and to reduce the number of cells; see Figure 2.

Figure 2: Cross section of the four proposed space discretization schemes for the domain (M1, M2, M3, and M4).

The nodalisations have the same number of axial levels (176). Geometrical aspects of the spatial discretization are summarized in Table 2.

3.1.4. Mesh Sensitivity Analysis. The case PSBT 1.6222 has been simulated using different meshes (M1 and M4). In Figure 3 the axial void fraction profile in the near wall region predicted by NEPTUNE-CFD for the two meshes is compared to each other. In Figure 4, the axial coolant temperature profiles predicted near the wall (T.near wall) and at the centre of the subchannel (T.center) for M1 and M4 are compared to each other. The largest difference can

Table 1: Test conditions for steady-state void measurements.

Run number Heat flux (W/cm2) T. inlet (K) Pressure (MPa) AT (K) Mass flow (kg/s)

1.2211 194.34 568.4 15.01 46.9 0.3248

1.2223 150.72 592.6 15.01 22.6 0.3248

1.2237 129.56 602.6 15.03 12.9 0.3248

1.4325 129.13 526.8 10.03 57.6 0.1487

1.4326 129.77 541.8 10.01 42.4 0.1487

1.6222 107.75 477.2 5.0 59.9 0.1487

Table 2: Discretization details of the four meshes applied in the study of the PSBT benchmark.

Subchannel 1 (M1) Subchannel 2 (M2) Subchannel 3 (M3) Subchannel 4 (M4)

Number of cells 182679 147429 112179 77345

First cell near the wall (mm) 0.1 0.2 0.3 0.4

Number of cross cells (X direction) 25 21 17 14

Axial height Z (m)

O VF near wall (M1) O VF near wall (M4)

Figure 3: Axial void fraction (VF) profile in the near heated wall region for case 1.6222.

Axial height Z (m)

O T.near wall (M1) - Txenter (M1)

O T.near wall (M4) —" T.saturation ---T.center (M4)

Figure 4: Axial water temperature profile in the near heated wall region and in the subchannel center for case 1.6222.

be observed for the void generation near the wall, where refined meshes generate more steam than coarse meshes. NEPTUNE-CFD predicts a fast axial void fraction increase using Ml (fine meshing) than using M4 (coarse meshing). For Ml, the maximal void fraction is calculated at 0.55 m, while for M4 this value shifts to 1 m elevation.

The heat flux partitioning at the wall is commonly divided into 3 fluxes according to Kurul and Podowski [13]. This model releases the heat flux from the heated wall into the water phase, and then it is subdivided into evaporation, quenching, and convective flux. Hence, the water phase receives all the heat fluxes. In the refined mesh (M1) this is problematic for the water, which cannot dissipate the

heat fast enough; hence some small peaks appear in the water temperature once the cell has been filled with vapour (Figure 4). This phenomenon is reproduced also in the coarse mesh but is not that severe. To avoid it, the wall heat flux is solved according to the four-flux model described in Section 2. The heat flux is redirected to the steam following a CHF criterion based on the void concentration on the near wall region. This measure can relax the water temperature, but if the heat flux is high enough, the problem is reproduced in the steam temperature. The run number 1.2211 has been selected to test the influence of the four different meshes on the NEPTUNE-CFD predictions. A constant bubble diameter of 0.1mm has been assumed for all mesh analysis. The

618 616 614 -

g 612 -

ÍU Ui

I 610 -

J 608 -

| 606 -604 602 600

CÀ 8 6

Ù> d ?

S 1 8

- \ Âf

0.0015

0.003 0.0045 X (m)

0.0075

O M1 A M2 M3

Figure 5: Water temperature in the near wall region at 1.4 m elevation as function of the mesh size for 0.1 mm bubble diameter.

location selected to visualize the local values of temperatures and void fraction (VF) is illustrated in Figure 1(b); this is a 8.32 mm distance between heated walls at Z = 1.4 m. The temperatures near the wall reach saturation for all cases as illustrated in Figure 5. A coarse nodalisation generates flatter profiles than finer meshes, especially regarding the void distribution calculation as shown in Figure 6. There, the steam production is higher for the smaller cells.

For the comparison of four different meshes, NEPTUNE-CFD predicts large gradients of void fraction in the near wall region for the refined meshes (M1) and (M2). These void fraction differences affect the calculated average void fraction over the cross-sectional area at the measurement position. The water temperature evolution calculated for the different mesh configurations at two axial locations can be observed in Figure 7. The first location is at centre of the subchannel, the second is in near wall heated region (first cell near the wall). At the centre of the subchannel, there are no big differences in the temperatures; only close to the outlet there are some deviations. In the near heated wall region locally the difference can reach 6 degrees K between refined mesh (M1) and coarse mesh (M4).

The liquid temperature profile is the combination of several phenomena. The liquid phase near the wall is heated by the wall heat flux, which is divided into convective, evaporation, and quenching heat flux. Once the bubbles are generated, they migrate and condense within subcooled liquid in the core of the flow and hence heat the liquid. The molecular and turbulent heat fluxes inside the liquid phase also modify the temperature profile.

0.55 -

0.00 0.000

o o

o a o a

A° °A

O O

A A

" + ° .......+

¿O OA

aO O;

<> O A + O <>

Ao O A

+ A ( O ■ ' 3 OA +

C M1 + M3

0.005 X (m)

A M2 O M4

Figure 6: Local void fraction (VF) at 1.4 m elevation as function of mesh size for 0.1 mm of bubble diameter.

620 615 610 605 600 595 590 585 580 575 570 565 -0.

03 0.17 0.37 0.57 0.77 0.97 1.17 1.37

Axial height Z (m)

— Ml_T_nearwall

■ - M2_T_near wall

■■■ M3_T_nearwall

-- M4_T_near wall

1 Ml_T_center

M2_T.center M3_T.center M4_T.center

Figure 7: Water temperature at the near wall region and at the centre of the subchannel of the case 1.2211 as function of mesh size.

The results computed by NEPTUNE-CFD are summarized in Table 3. The experimental value for the case 1.2211 is 3.8%. The first three nodalisations overpredict the averaged void fraction and the last one slightly underpredicts it. Nevertheless, this experimental value has an error of 3% in the VF measured and is not clear which discretization provides better results. The computational time for the coarse M4 mesh is around two hours; by applying more refined nodalisation, this time increases, reaching around seven hours ofcalculation for the mesh M1. Some of the cases selected to be simulated have rather high VF concentration at the measurement location; for example, case 1.4326 has a 53.1% VF. Case 1.6222 has a

Table 3: Results for the local and averaged VF for the different meshes applied for the simulation of case 1.2211.

Nearest heated wall cell VF (%)

Average VF (%) at measured cross section (Z = 1.4 m)

30.6% VF concentration, and as was illustrated in Figure 3, the first cell near the heated wall is quickly filled with steam. This amount of gas leads to overheating problems in the water and steam phase. The coarse mesh M4 can shift this negative effect in the water and steam temperatures to higher locations and in some cases avoid them.

The choice of the coarse mesh to perform the rest of the simulations has been made to preserve the numerical stability while solving the heat transfer problem. In addition, the coarse mesh (M4) is producing maximum y+ values around 300, which makes it still valid for the application of the selected k-e turbulence model.

Ml 67 8.1 s?

M2 57 6.7 (re 612

M3 40 4.6 tra 610

M4 15 3.1 S et

is 608 -

0.001 0.002 0.003

Distance from the wall X (m)

Diameter 0.1 mm Diameter 0.2 mm

Figure 8: Water temperature in the near heated wall region at 1.4 m elevation. Comparison for different bubble diameters. Case 1.2211.

3.1.5. Bubble Size Sensitivity Analysis. In the previous simulations, a constant bubble diameter (0.1 mm) has been applied. Incidence of other bubble diameters or the selection of an interfacial area equation (IAE) for the simulation is discussed in this subchapter. For the mesh M4 and the case 1.2211 previously studied, three different configurations for the IAC are considered. The first is by applying the previous 0.1mm constant diameter. In the second configuration the diameter is increased to 0.2 mm. The third applies one IEA, described by (10). The bubble diameter influences the simulation results. Therefore, the VF at the measurement location varies depending on this parameter and different results can be obtained. When the code especially is dealing with simulations where the main phenomena are subcooled boiling and low VF are expected.

For the three configurations proposed, the water temperature in the near wall region is illustrated in Figure 8. Here, rather good agreement between cases can be observed. The local bubble size calculated by the IAE is shown by Figure 9. The SMD calculated in this case is much lower, around 0.07 and 0.03 mm, compared with the other cases (0.1 and 0.2 mm). By calculating this small bubble size with an IAE, the IAC is higher compared to the other cases. As a consequence, the heat and mass transfer is higher and there are two important effects. First, the condensation into subcooled liquid is stronger, decreasing the VF within subcooled regions. Second, the boiling in the superheated near wall region is higher. These effects are illustrated in Figure 10. Here, the VF profile illustrates that steam bubbles are nucleated at the heated wall surface and condense in the subcooled liquid in the core of the flow. Larger bubbles (0.2 mm) do not condense so easily leading to a higher VF concentration in the bulk flow region. Concerning the impact of the bubble size in the steam velocity (Figure 11), slightly

O o r O

0 0.001 0.002 0.003 0.004

Distance from the wall X (m)

Figure 9: Bubble size distribution in case of using the IAE. Case 1.2211.

higher velocities are registered for the steam in the case of 0.2 mm bubble diameter in the bulk region. This can be explained as a consequence of a higher VF concentration in this region.

In Figure 12 the axial VF evolution in the centre of the subchannel is shown. In the centre of the subchannel the VF generated by the 0.2 mm bubbles is clearly higher due to the reasons previously explained. The void calculated by the IAE with bubblesof0.03 mmofdiameterislower duetothe strong condensation into subcooled liquid.

Table 4: Comparison of predicted and measured void fraction.

Run number

PSBT VF (%)

NEPTUNE-CDF VF (%)

Relative error (%)

ANSYS CFX 12.1 SST VF (%)

Relative error (%)

1.222300 1.223700 1.221100 1.432600 1.432500 1.622200

31.100 44.000 3.800 53.100 33.500 30.600

20.280 31.850 3.100 58.890 40.470 41.620

-34.790 -27.610 -18.420 10.900 20.810 36.010

21.200 29.100 11.700 50.600 34.300 26.600

31.833 33.864 -207.895 4.708 -2.388 13.072

í Í3 c

0.01 -

□ □

o o o □ c 8 <e □ -ê—

0.001 0.002 0.003

Distance from the wall X (m)

O Diameter 0.1 mm □ Diameter 0.2 mm

Figure 10: VF in the near heated wall region at 1.4 m elevation. Comparison for different bubble diameters. Case 1.2211.

0.001 0.002 0.003

Distance from the wall X (m)

O Diameter 0.1 mm □ Diameter 0.2 mm

Figure 11: Steam velocity in the near heated wall region at 1.4 m elevation. Comparison for different bubbles diameters. Case 1.2211.

The pressure drop in the channel is calculated also for the different bubble sizes. The simulated results are presented in Figure 13. Only the case predicted with the IAE exhibits higher values. In this case the concentration of bubbles in the near wall region is bigger leading to an increment of the pressure drop.

By using an IAE there are more parameters to control. One of the most important measures is to clip the value of the minimum bubble diameter. Very small bubble size can lead to numerical instabilities regarding the nondrag forces applied like the added mass force or the turbulent dispersion force. For the simulations the smallest bubble diameter of 0.01 mm is allowed in the computational domain.

3.1.6. Selected Results. The VF calculated by NEPTUNE-CFD for the six cases selected from the PSBT database is illustrated in Table 4 together with the experimental measurements. These simulations have been performed with a single IAE and the mesh M4.

Three experiments are overpredicted (1.2223,1.2237, and 1.2211) and the other three are underpredicted (1.4326,1.4325, and 1.6222) by NEPTUNE-CFD. In Table 4 the relative error between the experimental data and the computed results is shown. The simulations performed show a variation of the relative error from -34% to 36% compared with the measured PSBT data. It must be noted that the two-phase flow models of CFD codes are still under development, and hence the deviations of the predictions from the measured data are still considerable. In Table 4, the predicted void fraction by ANSYS CFX 12.1 using the shear stress turbulence model (SST) and the respective deviations from the experimental data are also given. These values were taken from [20] and they confirm the trends predicted by NEPTUNE-CFD that the deviations may be huge depending on the modelling being used to describe the turbulence effects among others.

The axial VF profile for each case is illustrated in Figure 14; the experimental measure and its measurement error (3%) are also included.

Axial height Z (m)

----- IAE

- Diameter 0.1 mm

---Diameter 0.2 mm

Figure 12: Axial VF profile at the centre of the subchannel. Comparison for 3 different bubbles diameters. Case 1.2211.

1.504E + 07

1.504E + 07 X:

\ : \ \

1.503E + 07 - V \

& 1.503E + 07 - \

I \ \ '

£ 1.502E + 07 - \

MH s ■ v

1.502E + 07 - !

.....\V.

1.501E + 07 -

.......

1.501E + 07 -I-i-i-i-i-i-i-

-0.03 0.22 0.47 0.72 0.97 1.22 1.47 Axial height Z (m)

----- IAE

Diameter 0.1 mm ---Diameter 0.2 mm

Figure 13: Pressure drop calculated for the 3 different bubble diameters. Case 1.2211.

3.2. The NUPEC BWR Full Size Fine

Mesh Bundle Test (BFBT) Benchmark

3.2.1. Scope and Description of the Benchmark. The BFBT void distribution benchmark [18] was made available by the Nuclear Power Engineering Corporation (NUPEC). It is one of the most valuable databases identified for thermal hydraulic modelling. The NUPEC database includes subchannel void fraction and critical power measurements in a representative (full-scale) BWR fuel assembly (FA). The high resolution and high quality of subchannel void fraction data encourage advancement in understanding and modelling of complex two-phase flow in real bundles and make BFBT experiments valuable for the NEPTUNE-CFD multiphase models validation. The benchmark consists of two parts: void distribution benchmark (Phase I) and critical power benchmark (Phase II). Each part has different exercises including simulations of steady-state and transient tests as well as uncertainty analysis. An exercise from phase I has been selected for the validation of NEPTUNE-CFD. The transient tests performed in the frame of this benchmark represent the thermal hydraulic conditions that may be encountered during a postulated BWR turbine trip transient without bypass and a recirculation pump trip. From this postulated transient scenarios important thermal hydraulic parameters are derived for the test, such as the evolution of the pressure, total bundle power, mass flow, and radial and axial power profile, which serves as initial and boundary conditions for the CFD simulations.

The test section of the experiment is a full sized 8 x 8 BWR FA (Figure 15), with sixty electrically heated rods

(12.3 mm diameter and 16.2 mm rod pitch) and one water channel (34 mm diameter). The electrically heated section is 3708 mm high; the heaters are surrounded by an insulator (boron nitride) and by the cladding (Inconel 600). An X-ray densitometer is used to measure the averaged void fraction at three different axial levels from the bottom from the heated section (Z = 0.67 m, Z = 1.72 m, and Z = 2.7 m). What the BFBT experimental data provides is the evolution in time of the averaged VF at the three axial levels mentioned.

To reproduce in the test bench the turbine trip without bypass and the recirculation pump trip, the BFBT database provides the evolution of the water mass flow rate, the system pressure, and the power. This data is used as boundary condition for NEPTUNE-CFD. In Figures 16,17, and 18 the evolution of the outlet pressure, the mass flow rate, and the power during the 60 seconds transient is given for both scenarios. The experiment was performed with a uniform axial power shape. The radial power shape is described by Figure 19. The water inlet temperature remains constant at 552°K during the experiment.

3.2.2. Applied NEPTUNE-CFD Models. The numerical simulation is performed with the models explained in Section 2; fc-e turbulence model for the liquid phase is used with the two-phase modified wall function. Local equilibrium turbulence model [2] is considered for the dispersed phase. The libraries for the liquid properties are provided by the IAPWS data [19]; furthermore steam close to saturation condition is assumed. For heat transfer through the gas/liquid interface, a thermal phase change model is applied. The heaters are modelled in terms of a heat flux boundary condition. At the inlet,

0 0.3 0.6 0.9 1.2 1.5 1.8 Z (m)

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0 0.3 0.6 0.9 1.2 1.5 1.8 Z (m)

—e— Sim. ■ Exp.

0 0.3 0.6 0.9 1.2 1.5 1.8 Z (m)

0 0.3 0.6 0.9 1.2 1.5 1.8 Z (m)

—0— Sim. ■ Exp.

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0 0.3 0.6 0.9 Z (m)

1.2 1.5 1.

0 0.3 0.6 0.9 Z (m)

—e— Sim. ■ Exp.

1.2 1.5 1.

Figure 14: Axial VF evolution for each case compared against experimental data.

a mass flow rate is set and the pressure is imposed at the outlet. The temporal evolution of mass flow rate, power, and pressure are taken from the experimental conditions (Figures 16,17, and 18).

The time step is adaptive depending on the Courant number. The time step width ranges from 1 to 3 ms. The drag and nondrag forces, lift, added mass, and turbulent dispersion force, are computed for the simulation.

3.2.3. Modelling of the Rod Bundle. The nodalisation of the PSBT test section was done following the best practice guidelines for the use of CFD codes. The nearest wall heated region cell has a constant width of 0.3 mm. The mesh is composed of 135 axial levels and 12 cross cells in each subchannel. Globally the NEPTUNE-CFD nodalisation has 211928 cells. The maximum y+ values are located close to the outlet. They oscillate between 300 and 400 depending on the transient conditions. These high values correspond to the high velocities locally achieved by the steam.

Taking into account the radial symmetry of the FA only 1/8 of the fluid domain is modelled with two symmetry planes (Figure 15). This decision is taken although being aware that the radial power distribution has no axial symmetry according to Figure 19, but larger models penalize the computational time. For the initialization of the simulations, the power is increased gradually from 0 to nominal values which are reached after 6 seconds. The simulations are conditioned by a large amount of void present in the domain. Numerical simulations find convergence problems in very refined meshes with boiling flow, especially with respect to the water and steam temperature close to the heated wall region. To set the heat flux defined by the experimental measurements (Figure 18) two different configurations are investigated. The first option is to locate the flux at the rod wall of the fluid domain using NEPTUNE-CFD as standalone, Figure 20(a). The second option is to model the clad and insulator of the test bench, where the heat flux is placed at the inner diameter of the insulator's surface, Figure 20(b). In both cases the power is set as superficial heat flux (W/m ).

Outlet

Top heated - section Z = 3.7 m

VF measured at axial level 3 S Z = 2.7 m

VF measured at axial level 2 Z = 1.72 m

VF measured at axial level 1 Z = 0.67 m

Bottom heated

section Z = 0 m / i ! 1 / 1 / r

i- 1/ j

¡La !

AY • Z

Cross section of the 8x8 FA

Water velocity (m/s)

Numerical model with the measure axial level colored

Outlet

Symmetry planes (d)

6.07838 fe 6

"E 4 3

1.43145

VF measure locations

= 2.7 m

,Z = 1.72 m

,Z = 0.67 m

Figure 15: (a) Test section with the location of the X-ray densitometers and the heated section. (b) Cross section of the tested rod bundle and the selected portion to be modelled. (c) Numerical model of the simulated domain. (d) Water velocities calculated by NEPTUNE_CFD.

NEPTUNE-CFD is not able to solve the thermal problem for the solid domain by itself. It has to be coupled with the heat conduction tool SYRTHES [21].

The averaged void fraction for the cross-sectional area of the FA is calculated each 0.02 seconds at the three axial levels (Z = 0.67 m, Z = 1.72 m, and Z = 2.7 m) from the bottom of the heated section. Figure 15 shows the location of the measured sections. The original experimental measurements were considered not to be accurate enough and a correction factor has been proposed to be applied; see [22]. D0 is the original experimental data in percentage of void, and D1 is the resulting data after the application of the correction factor. The new corrected data is obtained according to the following:

1 (-0.001 *D0 + 1.167)'

3.2.4. Selected Results for the Turbine Trip without Bypass Experiment. Computed void fraction results corresponding to two configurations are illustrated. The first is by modelling the problem with NEPTUNE-CFD standing alone and the

second configuration takes into account the wall thermal inertia effect by adding the SYRTHES code to NEPTUNE-CFD.

Figure 21 illustrates with different figures the local VF distribution calculated by NEPTUNE-CFD and SYRTHES in the domain for different time steps during the transient: at second 7, at second 11.5 (during the maximum void concentration), at second 20, at second 30, and at second 50, where the combination of power and mass flow decreases the void to the minimum values of the simulation. Furthermore the location of the measured axial levels is shown.

The comparison between the evolution of the calculated void averaged calculated by NEPTUNE-CFD and SYRTHES and the experimental data is shown in Figure 22 for the three different axial levels; the experimental measurement error of the mean averaged VF (±2%) is also included.

At this point, it is important to remark that the models implemented regarding interfacial interactions like mass, heat transfer, and momentum are designed for bubbly flow only. The water is set as continuous phase and the vapour is set as the disperse phase; this condition is valid for bubbly flow. When the amount of void in the fluid exceeds 60%,

- Pump trip

2G 30 4G

Time (s)

------- Turbine trip

Figure 16: Test section outlet pressure evolution measured for turbine trip and recirculation pump trip experiment.

7.G 6.5 6.G 5.5 £ 5.0Г 4.5

& 4.G-З.5 ЗЮ -2.5 2.G

Time (s)

- Pump trip

------- Turbine trip

Figure 18: Test bench measured power evolution for turbine trip and recirculation pump trip experiment.

S зo -

20 30 40

Time (s)

- Pump trip

------- Turbine trip

Figure 17: Test section mass flow rate evolution measured for turbine trip and recirculation pump trip experiment.

the regime is no longer bubbly flow. Therefore, the results of the second and third elevation cannot be considered reliable and more developments in the modelling are necessary to properly describe the flow at those elevations. The description of different flow regimes in the same domain requires the definition of a change of continuous and dispersed phase at each location. Nevertheless the code is able to trace the tendency of the void generation in line with the experimental data even in those flow regimes.

1.15 1.З0 1.15 1.З0 1.З0 1.15 1.З0 1.15

1.З0 0.45 0.89 0.89 0.89 0.45 1.15 1.З0

1.15 0.89 0.89 0.89 0.89 0.89 0.45 1.15

1.З0 0.89 0.89 0.89 0.89 1.15

1.З0 0.89 0.89 0.89 0.89 1.15

1.15 0.45 0.89 0.89 0.89 0.89 0.45 1.15

1.З0 1.15 0.45 0.89 0.89 0.45 1.15 1.З0

1.15 1.З0 1.15 1.15 1.15 1.15 1.З0 1.15

Figure 19: Normalized radial FA power coefficients for turbine trip and recirculation pump trip experiment in the context of the BFBT experiment.

Figure 23 illustrates the relative error of the computed void fraction against the measurements for the first axial level (Z = 0.67 m) together with the normalized values of the pressure mass flow rate and power of the transient. In Figure 23, it can be observed how the code struggles especially when the mass flow rate decreases. The code has a constant relative error around -10% during the rest of the transient. At the end of the transient, the predicted void fraction agrees very well with the data at all three axial elevations.

The coupling of NEPTUNE-CFD with SYRTHES leads to a significant improvement due to the accurate steam temperature calculation. In addition SYRTHES contributes to relax the temperature calculation and the solid-liquid interface increasing the time step and accelerating the simulation.

Figure 24 illustrates the evolution during the turbine trip simulation of the local temperatures of the water and steam together with the saturation and clad temperature.

Heat flux

(W/m2 ) Inner surface of rod 4

Fluid domain, numerical model

Fluid domain solved by NEPTUNEXFD

Solid domain solved by SYRTHES

Cladding, numerical model

Insulator, numerical model

Heat flux (W/m2 ) Inner surface of the insulator

Fluid domain, numerical model

(a) NEPTUNEXFD

(b) NEPTUNEXFD + SYRTHES

Figure 20: (a) NEPTUNE-CFD model; (b) NEPTUNE-CFD+SYRTHES model; green ring represents the insulator and the clad according to the real benchmark geometries.

t = 7 s (b)

t = 11.5 s (c)

t = 20 s (d)

t = 30 s (e)

0.500 n

Figure 21: (a) Location of the measured axial levels. Local VF distribution calculated by NEPTUNEXFD at different times for the turbine trip scenario: (b) second 7, (c) second 11.5, (d) second 20, (e) second 30, and (f) second 50.

The temperatures from Figure 24 are taken from one point near the wall region of the rod 4 at the end of the heated section, where the temperatures are expected to be the highest in the domain. By applying NEPTUNE-CFD and SYRTHES during the power peak, the steam temperature remains as maximum at 3 degrees oversaturation, while in other simulations without solving the thermal wall problem (NEPTUNE-CFD standalone), this temperature can reach locally hundreds of degrees oversaturation which is physically

not correct. On the other hand if the heat flux is defined at the insulator, the water and steam temperatures are calculated in a better way since SYRTHES is providing the wall temperature. Thus, the heat flux at the solid-fluid interface is calculated according to the following considering the wall

and fluid temperature (Tw, Ty) and a heat transfer coefficient:

85 80 75 70 65

£ 50 £ 45

Z = 2.7 m, experimental VF

Experimental error

35 30 25 20 15 10 5 0

7 11 15 19 23 27 31 35 39 43 47 51 55 59 Time (s)

Exp.3 3 ErrUp T1

3ErrDown

Figure 22: Comparison of the BFBT VF data and its measurement error (±2%) with the VF evolution predicted by NEP-TUNE/SYRTHES at three axial levels during the turbine trip experiment.

0.9 0.8 0.7 0.6 0.5 0.4 0.3

7 11 15 19 23 27 31 35 39 43 47 51 55 59

Time (s)

Pressure evolution Power evolution

Mass flow evolution Relative error

Figure 23: Normalized values of power, mass flow, and pressure evolution provided by the turbine trip experimental data together with relative error between VF predicted by NEPTUNE-CFD/SYRTHES and experimental VF data at axial level Z = 0.67 m.

where u* is the wall friction velocity and T+ is the nondimen-sional liquid temperature. Hence the heat flux from the wall to the steam phase is regulated avoiding excessive overheating. The wall temperature increases at those locations with high void fraction. In Figure 24, the local temperature jumps between the wall and the liquid oscillate from 20 to 30 degrees K during the transient, and when the power peak occurs, this difference rises up to 90 degrees K.

3.2.5. Selected Results for the Recirculation Pump Trip Experiment. Using the experience acquired from the turbine trip simulation, another exercise of the BFBT database is simulated. For this case, only the simulation combining NEPTUNE-CFD and SYRTHES is performed.

7 11 15 19 23 27 31 35 39 43 47 51 55 59

Time (s)

Water_temperature - Steam_temperature

Sat.temperature Clad temperature

Figure 24: Saturation, steam, clad, and liquid temperature calculated by NEPTUNE-CFD/SYRTHES. Comparison against the steam temperature calculated with NEPTUNE-CFD stand alone.

95 90 85 80 75 70 65 £ 60 £ 55 50 & 45 § 40 % 35 30 25 20 15 10 5

7 11 15 19 23 27 31 35 39 43 47 51 55 59 Time (s)

Figure 25: Comparison of the BFBT VF data and its measurement error (±2%) with the VF predicted by NEPTUNE/SYRTHES at three axial levels during the recirculation pump trip experiment.

The comparison between the experimental data provided by the BFBT database and the computed void fraction of NEPTUNE-CFD/SYRTHES is illustrated in Figure 25.

For this case the prediction fits well with the experimental data during the void increase between seconds 11 and 13. At this time there is no power peak and the void rises due to the decrease of mass flow. At this point the simulation is underestimating the experimental data only at the first axial level. For the lapse of time between seconds 15 and 42, an overprediction occurs for all axial levels. At the last third of the transient when the mass flow increases up to nominal conditions, the void predicted is slightly underestimated for axial levels 1 (Z = 0.67 m) and 2 (Z = 1.72 m), while for axial level 3 (Z = 2.7 m) a good agreement is achieved.

Figure 26 shows the local liquid temperature distribution for the 3 axial levels of the measurements at second 30 of the recirculation pump trip simulation. For the first axial level, the temperature is 3 degrees overheated in the near

562 °

- 561 Й

i 560 P 559

Figure 26: Local water temperature at 3 different axial levels, (a) 0.67 m, (c) 1.72 m, and (c) 2.7 m. Calculated by NEPTUNE-CFD/SYRTHES for the recirculation pump trip (second 30).

- 0.6 ja

Figure 27: Local void distribution at 3 different axial levels, (a) 0.67 m, (c) 1.72 m, and (c) 2.7 m. Calculated by NEPTUNE-CFD/SYRTHES for the recirculation pump trip (second 30).

Steam temperature at

the center of the subchannel with time

scale returning to

saturation 0.05 s

Steam temperature at

the center of the subchannel with time scale returning to saturation 0.01 s

Steam temperature in the bulk remains at saturation

579 577 575 g 573 1 571 3 569

£ 565 563 561 559

7 11 15 19 23 27 31 35 39 43 47 51 55 59 Time (s)

--Saturation_temperature

--Steam_bulk_ temperature

--Steam.walL temperature

--Steam_walL temperature _0.001

Figure 28: Local steam temperatures evolution at the bulk and near the heated wall calculated by NEPTUNE-CFD/SYRTHES during the recirculation pump trip. Comparison of two different time scales returning to saturation for the steam.

wall region and subcooled in the centre of the subchannels. In upper locations the steam concentration in the near wall region is above 80% (Figure 27), and all the heat flux is transferred into the steam relaxing the liquid temperature.

In Figure 28, the local temperature evolution is shown for two locations: the near wall region and the bulk of the fluid in the centre of one subchannel. Both locations are Z = 3.69 m high from the beginning of the heated section. The largest deviation from saturation takes place for the steam temperature in the near wall region. It can reach locally 17 K oversaturation. This steam temperature has different values by decreasing the time scale returning to saturation (8). This time scale should be lower than the time step to work properly. In this figure a time scale of 0.05 seconds is compared with a time scale of 1 ms. The temperature difference between those two settings is about 10 degrees. The lower time scale helps to maintain the steam temperature close to saturation in the near wall region. The temperature of the liquid in the near wall region is slightly oversaturated and it is not affected by the different time scales of the steam heat transfer coefficient. This parameter is not affecting significantly the steam generation or the temperature of the steam and the liquid at the bulk region which remains at saturation.

4. Summary and Conclusions

The validation basis of NEPTUNE-CFD has been extended by using LWR-relevant transient test data obtained in both a PWR-specific (NUPEC PSBT) and a BWR-specific (NUPEC BFBT) test facility. It is the first time that the two-phase

flow models of NEPTUNE-CFD have been validated using transient data obtained for tests representing the LWR plant conditions of postulated transients like a turbine trip or a recirculation pump trip. This validation process has clearly shown the status of the two-phase flow modelling of NEPTUNE-CFD. It is possible to summarize the global conclusions about the NEPTUNE-CFD capabilities, identifying its strengths and weaknesses.

According to the interfacial exchange terms, the main flow regime currently implemented in NEPTUNE-CFD is the bubbly flow, which is good enough to describe the flow with a certain void concentration (<60%). However, in BWR the amount of void generated is beyond the limits of the bubbly flow, and other flow regimes (slug flow or annular flow) must be taken into account for an accurate description of safety-relevant phenomena. How to model the transition between flow regimes describing a whole flow map is an open issue for CFD codes.

A study of the sensitivity to the bubble size has been performed applying two constant bubble sizes (0.1 and 0.2 mm) and a variable size with an IAE for the PSBT experiments. The heat and mass exchange is strongly influenced by the IAC. Bubbles may rise close to the heated wall and eventually depart from it and migrate into subcooled liquid where they condense. Smaller diameters produce more IAC and more condensation. Hence, for smaller selected bubble sizes, reduced void fraction will appear in the domain. The smaller bubble size generated by the IAE produces more interfaces, and thus the condensation is stronger and the void predicted is less than the other two bigger bubbles selected (0.1 and 0.2 mm).

By using an IAE it is assumed that there is a single bubble size per cell. If all bubbles have locally the same size, they will condense at the same speed and their diameter will decrease. If a multisize model is applied, the small bubbles will decrease and collapse rapidly leaving the bigger ones, which increase the mean bubble size. Therefore, the assumption of a single bubble size can lead to an underestimation of the bubble mean size in this case, affecting the void generation.

NEPTUNE-CFD alone cannot solve heat conduction in solid domains. For this reason it is coupled with SYRTHES which is in charge of calculating the temperature in solid domains. The capabilities of the NEPTUNE-CFD coupled with SYRTHES have been demonstrated by the prediction of thermal inertia at the walls. By applying SYRTHES, the simulations of the turbine trip are in better agreement with experimental data compared with the application of NEPTUNE-CFD standalone. In addition, this heat conduction solver helps to reasonably control steam and water temperatures in the near wall region by computing a proper heat transfer coefficient at the solid/fluid interface. The use of SYRTHES has a positive contribution to the NEPTUNE-CFD prediction capabilities. However, it is not parallelized and hence it penalizes the computational time.

Due to the different mesh distribution applied for each code (tetra volumes for SYRTHES and hexa volumes for NEPTUNE-CFD), a careful mesh design is required to match the nodes at the solid/liquid interface and have a good energy balance at the interface.

The classical 3 fluxes formulation for the wall/fluid heat transfer assumes three fluxes: quenching, convection, and evaporation. If boiling occurs in this region and the amount of liquid decreases, the water enthalpy rises and the temperature can be several degrees above saturation. To avoid this problem, the four-flux model decomposition (see (25)) transfers all the heat flux into the steam when the void is above the 80% in the near wall region (condition frequently present near the wall). With this measure, the liquid temperature remains close to saturation, but now the steam can be locally overheated in the near wall region if the heat flux increases. The steam temperature problem can be fixed by increasing the near wall cell size. By this way, the y + values increase; it is important to maintain those values in a range according to the turbulence model selected. The turbulence model selected in this case is the two-equation based fc-e. The wall-adjacent cells must belong to the log region of the wall boundary layer (30 < y+ < 300) for this model to work properly. Other more sophisticated turbulent models have been tested, for instance, a 7-equation model, but no convergence has been obtained, mainly because they require a refined near wall region discretization to reach y + values below 10.

The near wall cells size results as an agreement between the thermal and momentum modeling. This agreement consists in solving the heat transfer problem properly without penalizing excessively the y+ values.

By setting a very refined axial nodalisation (more than 200 axial levels for the presented domains), the code struggles calculating the pressure field. It is recommended to enlarge the cells in the axial direction in the locations where large amounts of steam are generated.

According to the results obtained during the validation process, the subcooled and saturated boiling description can be considered sufficient enough in its range of applicability. Even if further developments are required, NEPTUNE-CFD has demonstrated tobea valid TH tool for the two-phase flow modeling in LWR applications.

List of Symbols

Hki: K:

Ja: kk: L:

Interfacial area concentration (m2/m3)

Wall thermal diffusivity (m2/s)

Liquid thermal diffusivity (m2/s)

Thermal capacity of the liquid (J/(mol-K))

Gas heat capacity at constant pressure (J/(mol-K))

Sauter mean bubble diameter (m)

Bubble detachment diameter (m)

Bubble detachment frequency (Hz)

Interfacial-averaged enthalpies (kJ)

Heat transfer coefficient

Jakob number

Thermal conductivity (W/(m-K)) Latent heat of vaporization (kJ/kg) Interfacial transfer of momentum Active nucleation sites density Nusselt number Liquid Prandtl number Peclet number

qc: Convective heat flux (W/m2) qe: Evaporation heat flux (W/m2) Quenching heat flux (W/m2) q^'k: Volumetric heat flux (W/m2) qki: Interfacial heat flux density Refc: Bubble Reynolds number St: Stanton number Tk: Temperature of phase k (K) T+: Nondimensional liquid temperature tq: Quenching time (s) u *: Wall friction velo city (m/s) Vk: Phase k velocity (m/s) WC: Wall-adjacent cell y+: y plus value a: Void fraction r^: Mass transfer condition (Kg) vk: Kinematic viscosity (m2/s) r: Time scale returning to saturation (s) 0nuc: Bubble nucleation source terms 0coa: Bubble coalescence source term 0brk: Bubble break-up source term As: Wall conductivity (W/(m-K)) pk: Phase k density (kg/m2).

Abbreviations

BFBT: Boiling water reactor full bundle test

BWR: Boiling water reactor

CEA: Commissariat de l'Energie Atomique

CFD: Computational fluid dynamics

CHF: Critical heat flux

DNB: Departure from nucleated boiling

EDF: Electricité de France

IAE: Interfacial area equation

IAC: Interfacial area concentration

IRSN: Institut de Radioprotection et de Sureté

Nucléaire

IAPWS: The International Association for the

Properties of Water and Steam

KIT: Karlsruhe Institute of Technology

NUPEC: The Nuclear Power Engineering

Corporation

RANS: Reynolds Averaged Navier Stokes

SMD: Sauter mean bubble diameter

SST: Shear stress turbulence model

OECD: Organisation for Economic Cooperation

and Development

PSBT: Pressurized water reactor subchannel and

bundle tests

PWR: Pressurized water reactor

VF: Void fraction.

Conflict of Interests

The authors hereby declare that no conflict of interests is present between them and the commercial entities mentioned in the context of the paper.

Acknowledgments

The authors thank the Program "Nuclear Safety Research" of

KIT for the financial support of the Research Topic "Multi-

physics Methodologies for Reactor Dynamics and Safety" and

the EU Project NURISP.

References

[1] A. Guelfi, D. Bestion, M. Boucker et al., "NEPTUNE: a new software platform for advanced nuclear thermal hydraulics," Nuclear Science and Engineering, vol. 156, no. 3, pp. 281-324, 2007.

[2] J. Lavieville, E. Quemarais, and M. A. M. L. Boucker, "Neptune-CFD user guide," 2010.

[3] J. Lavieville, E. Quemarais, M. Boucker, and S. Mimouni, Neptune-CFD V1.0 Theory Manual, 2006.

[4] M. Ishii and N. Zuber, "Drag coefficient and relative velocity in bubbly, dropet or particulate flows," AIChE Journal, vol. 25, no. 5, pp. 843-855,1979.

[5] N. Zuber, "On the dispersed two-phase flow in the laminar flow regime," Chemical Engineering Science, vol. 19, no. 11, pp. 897917,1964.

[6] M. Ishii, "Two-fluid model for two-phase flow," Multiphase Science and Technology, vol. 5, no. 1-4, pp. 1-63,1990.

[7] T. R. Auton, "The lift force on a spherical body in a rotational flow," Journal of Fluid Mechanics, vol. 183, pp. 199-218,1987.

[8] A. Tomiyama, H. Tamai, I. Zun, and S. Hosokawa, "Transverse migration of single bubbles in simple shear flows," Chemical Engineering Science, vol. 57, no. 11, pp. 1849-1858, 2002.

[9] M. Lance and M. Lopez de Bertodano, "Phase distribution phenomena and wall effects in bubbly two-phase flows," Multiphase Science and Technology, vol. 8, no. 1-4, pp. 69-123,1994.

[10] W. Yao and C. Morel, "Volumetric interfacial area prediction in upward bubbly two-phase flow," International Journal of Heat and Mass Transfer, vol. 47, no. 2, pp. 307-328, 2004.

[11] Q. Wu, S. Kim, M. Ishii, and S. G. Beus, "One-group interfacial area transport in vertical bubbly flow," International Journal of Heat and Mass Transfer, vol. 41, no. 8-9, pp. 1103-1112,1998.

[12] Y. Hsu, "On the size range of acting nucleation cavities on a heating surface," Journal of Heat Transfer, vol. 3, no. 84, pp. 207216, 1962.

[13] N. Kurul and M. Z. Podowski, "Multi dimensional effects in forced convection subcooled boiling," in Proceedings of the 9th International Heat Transfer Conference, p. 21, Jerulasem, Israel, 1990.

[14] H. C. Unal, "Maximum bubble diameter, maximum bubble-growth time and bubble-growth rate during the subcooled nucleate flow boiling of water up to 17.7 MN/m2," International Journal of Heat and Mass Transfer, vol. 19, no. 6, pp. 643-649, 1976.

[15] Y. Egorov and F. Menter, "Experimental Implementation of the RPI wall boiling model CFX-5.6," Technical Report ANSYS/TR-04-10, ANSYS GmbH, 2004.

[16] J. Weisman and B. S. Pei, "Prediction of critical heat flux in flow boiling at low qualities," International Journal of Heat and Mass Transfer, vol. 26, no. 10, pp. 1463-1477,1983.

[17] A. Rubin, A. Schoedel, M. Avramova et al., "OECD/NRC Benchmark based on NUPEC PWR subchannel and bundle test (PSBT)," NEA-1849 ZZ-PSBT, 2010.

[18] D. Neykov, F. Aydogan, L. Hochreiter et al., "NUPEC BWR full-size fine-mesh bundle test (BFBT) Benchmark Volume I: specifications," Nuclear Science NEA/NSC/DOC, vol. 92, no. 64, 2005.

[19] Anon, Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, The International Association for the Properties of Water and Steam, Erlangen, Germany, 1997.

[20] E. Krepper and R. Rzehak, "CFD analysis of a void distribution benchmark of the NUPEC PSBT tests: Model calibration and influence of turbulence modelling," Science and Technology of Nuclear Installations, vol. 2012, Article ID 939561,10 pages, 2012.

[21] C. Péniguel and A. I. Rupp, SYRTHES 3.4—Manuel Theorique, ED F R&D, Paris, France, 2004.

[22] M. Gluck, "Validation of the sub-channel code F-COBRA-TF. Part II. Recalculation of void measurements," Nuclear Engineering and Design, vol. 238, no. 9, pp. 2317-2327, 2008.

Copyright of Science & Technology of Nuclear Installations is the property of Hindawi Publishing Corporation and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.