Multidimensional analysis of fission gas transport following fuel element failure in sodium fast reactor

Nuclear Engineering and Design 247 (2012) 136–146

Contents lists available at SciVerse ScienceDirect

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

Multidimensional analysis of fission gas transport following fuel element failure in sodium fast reactor Igor A. Bolotnov a , Steven P. Antal b , Kenneth E. Jansen c , Michael Z. Podowski b,∗ a b c

North Carolina State University, Raleigh, NC, United States Rensselaer Polytechnic Institute, Troy, NY, United States University of Colorado, Boulder, CO, United States

a r t i c l e

i n f o

Article history: Received 19 August 2011 Received in revised form 3 February 2012 Accepted 6 February 2012

a b s t r a c t Significant progress in several areas will have to be made to achieve the required technological and safety standards for future Gen. IV reactors, including both novel experimental methods (starting with separate-effect, then followed by integral experiments) and high performance computational models characterized by the necessary level of modeling detail and high accuracy of predictions. Furthermore, it is important that the experimental and theoretical/computational research complement each other, so that the results of measurements could be directly used for model validation purposes, whereas the results of simulations should provide input to identify modeling uncertainties and provide guidelines for prioritizing future experiments. The purpose of this paper is to present the modeling concept for mechanistic computer simulations of the injection of a jet of gaseous fission products into a partially blocked SFR coolant channel following localized cladding overheat and breach. A three-dimensional model of gas/liquid-sodium interaction has been developed based on a multifield modeling framework implemented in the NPHASE-CMFD code. The boundary conditions used as input to NPHASE-CMFD have been obtained by averaging the results of direct numerical simulations (DNS) performed using the PHASTA code. The novel aspects of the results discussed in the paper include the demonstration of advantages of using a multiscale approach to model local phenomena governing gas/liquid-sodium two-phase flow inside reactor coolant channels following cladding breach, as well as the observations about areas where future experiments are needed to improve the predictive capabilities of two-phase flow models. © 2012 Elsevier B.V. All rights reserved.

1. Introduction The overall objective of the study has been to develop an advanced multiscale approach to model local phenomena inside reactor fuel material, gas gap, cladding and coolant channel in a Gen. IV Sodium Fast Reactor (SFR), anticipated to occur during an accident caused by local coolant flow reduction. Since a general analysis of reactor accidents encompasses several phenomena of different nature, such as fluid mechanics and heat transfer, solid mechanics, fluid-structure interactions, fission product chemistry, a consistent method aimed at improving our understanding of mechanisms governing accident progression is via analyzing the individual phenomena in a parametric fashion, both in a standalone mode and including their coupling to other phenomena in

∗ Corresponding author. Tel.: +1 518 276 4000; fax: +1 518 276 2623. E-mail addresses: igor [email protected] (I.A. Bolotnov), [email protected] (S.P. Antal), [email protected] (K.E. Jansen), [email protected] (M.Z. Podowski). 0029-5493/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2012.02.004

the form of boundary conditions and/or interfacial interactions. As an example, an analysis reported by Bolotnov et al. (2010) can be mentioned here on the effect of fuel element heatup on cladding failure mechanisms and breach growth. The focus of the current work has been on studying multiscale aspects of fission gas release from failed fuel elements and its transport along reactor coolant channels for conditions corresponding to local channel blockage. The concept is based on using multiple computer codes to perform multiphysics simulations governing accident progression. The specific objective of this paper is to illustrate the proposed approach, and to show the advantages of using two different computer codes, NPHASE-CMFD and PHASTA, in a study aimed at three-dimensional (3-D) simulations of fission gas jet injection into a reactor coolant channel filled with liquid sodium. Channel blockage accidents can be caused by a local flow blockage by foreign objects or by mechanical failure/bowing of the hottest fuel pin. The main events characterizing accident progression include the following: • imbalance between power generation and heat removal rates,

I.A. Bolotnov et al. / Nuclear Engineering and Design 247 (2012) 136–146

137

Fig. 1. SFR fuel configuration used in present analysis: (a) hexagonal fuel assembly, (b) hottest fuel rod and the surrounding coolant channels, and (c), fuel pin (fresh fuel).

• fuel temperature increase, • fission product gas pressure increase, both inside and outside (in the expansion region) the fuel material, • local heatup, deformation and failure of stainless steel (SS) cladding wall, • injection of multicomponent fission product gas, possibly mixed with small fuel particles, into the coolant channel, • transport of sodium/fission-product mixture along coolant channels toward core outlet, • reactor power change due to fuel relocation (local reactivity effects).

The geometry of a typical hexagonal fuel assembly of SFR is shown in Fig. 1(a). The current analysis and simulations are focused on a single (hottest) fuel pin and the surrounding coolant channel. This is shown in Fig. 1(b). The initial phase of the accident under consideration is the heatup of the intact hottest cylindrical fuel element shown in Fig. 1(c). Whereas two different fuel materials have been considered for future SFRs, i.e., either oxidic or metallic, that latter has been used in the present due to its lower melting temperature and the impact on cladding failure and fission gas release. As a result of fuel heatup, cracks inside the fuel pin gradually develop and gaseous fission products get released into the gap between the fuel and the cladding at an accelerated manner. When the fuel material locally reaches the fuel melting temperature, the molten zone propagates outwards, thus accelerating the heatup of cladding. Eventually, a breach of the cladding wall occurs, and a jet of pressurized gas is injected into the coolant channel filled with molten sodium. As it is illustrated in Fig. 2, the overall accident propagation scenario involves several interacting phenomena which must be analyzed at different scales. The specific objective of this paper is to develop a mechanistic three-dimensional (3-D) model of the interacton between the jet of fission product gas and the liquid sodium flowing along the hot reactor channel. As indicated before, the main objective of the present work has been to study the fluid mechanics of phenomena governing fission gas injection into liquid sodium coolant at low-flow conditions. Consequently, parameters such as breach geometry, gas composition and injection rate, and others, have been used as input to the models. This, in turn, allowed one to study the effect of turbulence

Fig. 2. A schematic of fuel degradation and transport in SFR hot channel during channel blockage accident.

and fuel channel geometry on gas/liquid interaction during the postulated local core blockage accident. 2. Problem formulation For the purpose of the analysis, a section of the SFR fuel assembly shown in Fig. 1(b) has been selected, as shown in Fig. 3(a). A crosssectional view across the opening in cladding wall after the failure is shown in Fig. 3(b). It has been assumed for the purpose of the Reynolds-Averaged Navier Stokes (RANS) simulations that the jet

138

I.A. Bolotnov et al. / Nuclear Engineering and Design 247 (2012) 136–146

Fig. 3. Fuel rod geometry for side jet-injection problem: (a) top view and (b) side view.

Fig. 4. Three-dimensional (3-D) grid used in computer simulations of gas jet injection into liquid sodium.

(LES), detached eddy simulation (DES), direct numerical simulation (DNS)) flows. PHASTA (and its predecessor, ENSA) was the first unstructured grid LES code (Jansen, 1993) and has been applied to turbulent flows ranging from validation benchmarks (channel flow, decay of isotropic turbulence) to complex flows (airfoils at maximum lift, flow over a cavity, near lip jet engine flows and fintube heat exchangers). The PHASTA code uses advanced anisotropic adaptive algorithms (Sahni et al., 2006) and the most advanced LES/DES models (Tejada-Martinez and Jansen, 2005). Note that DES, LES, and DNS are computationally intensive even for single phase flows. This capability has been recently (Nagrath et al., 2006; Bolotnov et al., 2011) extended to two phase flows where we use the level set method to track the boundary between two immiscible fluids (either compressible – where we captured new instabilities in sonoluminescence, or incompressible – to study bubble coalescence and two-phase turbulence). PHASTA uses anisotropically adapted unstructured grids and its highly scalable performance on massively parallel computers has already been demonstrated (the code has shown perfect scaling out to 288*1024 IBM Blue Gene processors). The spatial and temporal discretization of the Incompressible Navier-Stokes equations within PHASTA has been described in Whiting and Jansen (2001). The strong form of the INS equations is given by Mass conservation

is symmetric relative to vertical plane crossing the direction of the jet inflow. This is shown on Fig. 3(a). Whereas the current geometry was deemed appropriate for the present analysis, aimed at establishing an accurate modeling and computational framework, its extension to account for a multiple fuel rod structure is planned for future work. The liquid sodium parameters used as a reference correspond to an elevated bulk temperature of 600 ◦ C (above the normal range of 500–550 ◦ C, but below the sodium melting temperature of 883 ◦ C) and a reduced velocity of 1.0 m/s. The fission gas injected through the side jet was assumed to have the properties of nitrogen at 700 ◦ C and 1.3 MPa. It is interesting to mention that the sodium flow conditions correspond to an average Reynolds number of about 35,000. A complete computational 3-D computational domain has been used. The corresponding unstructured grid, including the fission gas jet injection zone, is shown in Fig. 4.

ui,j = 0

3. Modeling and computational approach in the PHASTA code

ˆ ij = 2Sij = (ui,j + uj,i )

PHASTA is a parallel, hierarchic, high-order accurate (from the 2nd to the 5th order accuracy, depending on function choice), adaptive, stabilized (finite element) transient analysis flow solver (both incompressible and compressible), that has been developed at RPI. This approach has been shown by Jansen (1999) and Whiting and Jansen (2001) to be an effective tool for bridging a broad range of length scales in turbulent (RANS, large-eddy simulation

Using the Continuum Surface Tension (CST) model of Brackbill et al. (1992), the surface tension force is computed as a local interfacial force density, which is included in fˆi . The stabilized finite element formulation used in the PHASTA code (Whiting and Jansen, 2001) is based on the work of Taylor et al. (1998), Franca and Frey (1992). The predictive capabilities of PHASTA for turbulent flows have been extensively investigated before. As it was shown by Trofimova

(1)

Momentum conservation ui,t + uj ui,j = −ˆp,i + ˆ ij,j + fˆi

(2)

where  is density, ui is the ith component of velocity, pˆ is pressure, ˆ ij is the stress tensor, and fˆi represents body forces along the ith coordinate. For the incompressible flow of a Newtonian fluid, the stress tensor is related to the strain rate tensor, Sij , as (3)

I.A. Bolotnov et al. / Nuclear Engineering and Design 247 (2012) 136–146

et al. (2009) an excellent agreement has been observed between the results of PHASTA and the widely accepted law of the wall consistent with the experiments data of Laufer (1954). More recently, the results of PHASTA simulations for turbulent channel flows have been compared by Bolotnov et al. (2011) against the seminal results of Moser et al. (1999), and both the consistency and accuracy of the numerical predictions has been demonstrated. Based on the past applications of PHASTA, including but not limited to those referenced above, one concludes that for a wide range of turbulent flow conditions the results of this code’s simulations can be used as “virtual experimental data” of reference for CFD RANS-averaged models.

∂[(1 − ˛)c v¯ c ] + ∇ · [(1 − ˛)c v¯ c v¯ c ] ∂t = −(1 − ˛)∇ pc + (pc − pic )∇ ˛ + (1 − ˛)∇ · - c − (c −  ic ) · -

∇ ˛ + (1 − ˛)c g + Mic −  vic

(7)

∂ (˛d v¯ d ) + ∇ · (˛d v¯ d v¯ d ) ∂t

∇ ˛ + ˛d g + Mid +  vid

The NPHASE-CMFD code (Antal et al., 2000) is a robust computational multiphase fluid dynamics (CMFD) finite volume flow solver. The technology used by the NPHASE-CMFD code is the multifield model of multiphase flows. Separate equations for the conservation of mass, momentum and energy for each fluid/field have been developed as the solver’s framework. The governing equations are then ensemble-averaged, which allows the NPHASE-CMFD code to predict local non-steady-state multidimensional flow and heat transfer phenomena in two- or multiphase/multicomponent fluids. Key features of NPHASE-CMFD code include: • Finite-volume based solver • Use of unstructured grids with arbitrary element types • Capability to model an arbitrary number of fields (fluid components and/or phases) • Built-in mechanistic modeling, integrated with numerics • Improved robustness and numerical convergence • Modeling of gas/liquid interfaces using Level-Set method A typical form of conservation equations for mass, momentum and energy, respectively, of interpenetrating fluids can be written as (Podowski, 2009): (4)

∂(˛k k v¯ k ) t i + ∇ · (˛k k v¯ k v¯ k ) = −∇ (˛k p¯ k )+∇ · (˛k ␶ -¯ k ) + ˛k k g + Mk ∂t (5)

∂(˛k k e¯ k ) + ∇ · (˛k k v¯ k e¯ k ) ∂t t = −∇ · (˛k q¯ k ) + ∇ · [˛k (−p¯ k I + ¯ tk ) · v¯ k ] + ˛k k g · v¯ k + Eik - -

conservation equations for the liquid (continuous field, ‘c’) and gas (dispersed field, ‘d’) can be written as

= −˛∇ pd − (pd − pid )∇ ˛ + ˛∇ · - d + (- d − - id ) ·

4. Modeling and computational approach in the NPHASE-CMFD code

∂(˛k k ) + ∇ · (˛k k v¯ k ) = k ∂t

139

(6)

where ˛k is the volume fraction of field-k, k is the density of field-k, v¯ k is the velocity of field-k,  k is the mass transfer rate to field-k, p¯ k  t Re is the pressure in field-k, ␶ -¯ k = ␶ -¯ k + ␶ -¯ k is the total combined shear and turbulent shear stress for field-k, g is the gravity vector, e¯ k is the t k Re energy of field-k, q¯ k = q¯ k + q¯ k is the total heat flux for field-k, Mik is the interfacial force per unit volume exerted on field-k by the

other fields, and Eik is the corresponding interfacial energy transfer rate. Whereas the interfacial mass transfer is directly related to the net heat transfer rate at the interface, the interfacial transfer of momentum and energy also involve the effects of interfacial pressure and shear stress. In the case of gas/liquid bubbly flows, a two-field version of the above given model is used. In particular, the momentum

(8)

where ˛ = ˛d is the local void fraction and  =  d = −  c is the local volumetric evaporation rate. The model given by Eqs. (7) and (8) uses the Eulerian frame of reference and is based on the assumption that the individual fields can be treated as ‘equal partners’. However, it has been shown by Podowski (2009) that dispersed-gas/liquid flows do not follow, rigorously speaking, the interpenetrating-media concept. Fortunately, it has also been shown by Podowski (2009) that Eqs. (7) and (8) can still be used for gas/liquid flows, provided the following conditions are satisfied pd = pid = pic = pc = p

(9)

 = i = i =  = 

(10)

d

d

c

c

Multifield ensemble-averaged models normally require several closure laws for the various interfacial transport phenomena. In particular, the interfacial momentum transfer can be determined by expressing the overall interfacial force as a superposition of the individual forces governing specific transport mechanisms between the fields VM L TD Mid−c = −Mic−d = MD d−c + Md−c + Md−c + Md−c + · · ·

(11)

D VM VM where MD d−c = −Mc−d is the drag force, Md−c = −Mc−d is the vir-

TD tual mass force, MLd−c = −MLc−d is the lift force, MTD d−c = −Mc−d is the turbulent dispersion force. Other forces can also be included, such as the wall (lubrication) force. The analytic expressions for the individual forces have already been extensively documented, and can be found in Anglart et al. (1997), Drew and Passman (1999), Antal et al. (2005), and Podowski (2009). Another important modeling issue is concerned with two-phase flow turbulence. Typically, turbulence in multifield models used by CFD codes (such as Fluent, CFX, and others) is modeled using various versions of the k–ε model applied to the continuous field (Bertodano et al., 1990; Anglart et al., 1997; Antal et al., 2000, 2005; Antal and Podowski, 2003; Tiwari et al., 2006, 2009). In the case of the commonly used k–ε model(s), the conservation equations for turbulent kinetic energy and energy dissipation can be written as



˛c

Dk = ˛c ∇ · Dt

˛c

Dεc = ˛c ∇ · Dt

cT ∇k k





cT ∇ εc ε

+ ˛c (Pc − εc )

 + ˛c

εc (Pc C1˛ − C2˛ εc ) kc

(12)

(13)

where k is the turbulent kinetic energy, ε is the turbulent dissipation rate, cT is the turbulent viscosity of the continuous (liquid) field, Pc is the rate of turbulent energy production. The specific form of the expression for Pc depends on the particular version of the k–ε model, the High-Reynolds Number (High-Re)

140

I.A. Bolotnov et al. / Nuclear Engineering and Design 247 (2012) 136–146

Fig. 5. Stream-wise instantaneous velocity distribution for low Reynolds number turbulent flow (Re␶ = 178).

model (Jones and Launder, 1972) or the Low-Reynolds Number (Low-Re) model (Chien, 1982). The coupling between the turbulence model and the continuous-field momentum conservation equation in single-phase flows is introduced via the continuousfield turbulent viscosity term, given by t = C f

kc2 εc

(14)

where f = 1 for the High-Re model and f = 1 − exp(− 0.0115y+ ) for the Low-Re model (Chien, 1982). Since both versions of the k–ε model, as well as any other models based on using Reynolds stress averaging, are associated with substantial uncertainties, the model given by Eqs. (12) and (13) has been validated against a DNS solution. This is particularly important for situations where the turbulence distribution at the conduit inlet is used as a boundary condition in complex flow situations. Details in this regard are discussed in the next section. In the case of two-phase flows, the liquid-field turbulence may be augmented by the bubble-induced turbulence. The bubbleinduced kinematic viscosity of the liquid phase can be determined from the model proposed by Sato and Sekoguchi (1975) 2

c

=

Nd 





Cb db,k ˛k v¯ rel,k 

(15)

k=1

where Nd is the number of dispersed fields, Cb is the model coef  ficient, db,k is the average bubble diameter in the field-k, v¯ rel,k  is the magnitude of the relative velocity of the field-k. The Sato and Sekoguchi model has been successfully used before (Bertodano et al., 1990), but is important to recognize its limitations (in particular, as applied to non-bubbly flows). Past applications and validations of the NPHASE-CMFD code extend over a broad range of both single-phase and multiphase problems. Concerning the former, this code has been used to analyze flow and heat transfer in Gen. IV Supercritical Water Reactors (SCWR) (Gallaway et al., 2008) and in Very High Temperature GasCooled Reactors (VHTGR) (Gallaway et al., 2007). Multiphase flow applications include a study of steam injection into a large liquid pool (Nagrath et al., 2001), particulate flows in complex geometries (Tiwari et al., 2006, 2009), local effects (Antal et al., 2005) and oscillations (Antal and Podowski, 2003) in two-phase flow in advanced boiling water reactors, and bubble shape evolution in contact with solid walls (Wierzbicki et al., 2007). Recent studies on the validation of the NPHASE-CMFD code and the associated two-phase flow models include: a comparison against TOPFLOW data for multiplebubble-size gas/liquid flows in vertical conduits (Tselischeva et al., 2010a), predictions of void distribution in a horizontal pipe with 90 deg bend (Tselischeva et al., 2010b), and comparisons against the

data for void fraction and velocity distributions in gas/liquid flows (Shaver et al., 2011). 5. Validation of turbulence model It turns out that the highest level of uncertainty in the modeling of two- and multiphase flows is associated with the model of turbulence. Since in gas/liquid flows small-scale turbulence is practically limited to the continuous liquid field, a single-phase flow analysis based on the results of DNS simulations as a frame of reference, can be used to validate the RANS k–ε model(s). The purpose of this section is to demonstrate the usefulness of such an approach and to confirm the consistency and accuracy of the RANS turbulence models implemented in NPHASE-CMFD. In the example shown below, turbulence data has been transferred between the PHASTA and NPHASE-CMFD codes using a fully-developed solution of turbulent channel flow (Trofimova et al., 2009). The channel geometry and flow conditions were consistent with those of Moser et al. (1999). The channel flow was represented by rectangular non-dimensional geometry with the following parameters: • Mesh nodalization: N1 = 295; N2 = 129; N3 = 129. • Flow domain size: L1 = 4␲ (stream-wise); L2 = 2.0 (channel width); L3 = 4␲/3 (span-wise). • Re = u Ly /2 = 178.12, where u = (w /)1/2 is the friction velocity and  is the fluid kinematic viscosity. The flow was assumed to be periodic at the L1 and L3 boundaries, and there was no slip on the walls of the flow channel (i.e., the upper and lower L2 boundaries). Fig. 5 shows the stream-wise velocity distribution for a singlephase turbulent flow. As can be seen, PHASTA-DNS is capable of accurately determining the fluctuating velocity components of the solution at any point of the domain. A standard High-Re k–ε model discussed in Section 3 has been used in the NPHASE code for the first test case. In order to properly couple the PHASTA-DNS solution to the NPHASE-RANS problem, a grid for the same geometry has been created for NPHASE. This grid (the right side of Fig. 6) is quite coarse since the High-Re k–ε model (Jones and Launder, 1972) does not require a fine mesh and it also requires that the following conditions be satisfied simultaneously (Pope, 2000): y+ > 30, and y/ı < 0.3. The coarse NPHASE grid has been used to extract the coordinates of the cell centers in the inflow plane for use as monitor points in the PHASTA code. For the case shown on Fig. 6, there were only 6 points since the NPHASE grid was de facto two-dimensional. These

I.A. Bolotnov et al. / Nuclear Engineering and Design 247 (2012) 136–146

141

Fig. 6. PHASTA turbulent channel flow domain (left) outflow data is time-averaged to obtain inflow boundary conditions for NPHASE (High Reynolds number k–ε model case) domain (right).

20 18 16 14 12

U+

10 8 6 4 2 0

Fig. 7. Velocity components evolution at a single point of PHASTA-DNS (1 is streamwise velocity (u), 2 is normal to the wall (v) and 3 is the span-wise velocity (w) component).

monitor points have been supplied to the PHASTA code to record the velocity evolution at each of them. Fig. 7 shows the dependence of the velocity components at one of those points vs. the PHASTA time step number. The data obtained from the PHASTA run was post-processed to obtain the mean velocity, turbulent kinetic energy and turbulent dissipation rate using the following formulae: 1 ui (tj ) N N

Ui =

(16)

j=1

1  (ui (tj ) − Ui )2 N 3

N

k=

(17)

j=1 i=1

  1  ∂ui (tj ) ∂ui (tj ) N ∂xk ∂xk N

ε=

3

3

(18)

j=1 i=1 k=1

where ui (tj ) is the velocity component at a given point at time tj , Ui is the mean velocity component, ui ≡ ui − Ui is the fluctuating velocity component, N is number of time steps in consideration, k is the turbulent kinetic energy, ε is the turbulent dissipation rate, and  is the fluid kinematic viscosity. The mean velocity profile obtained by averaging the singlephase turbulent channel flow DNS performed by PHASTA agrees with law of the wall (Fig. 8) for y+ > 30: U+ =

1 ln y+ + C + 0

(19)

+ where U+ = U1 /u is the  dimensionless velocity, y = yu / is the w /1 is the friction velocity based on the wall coordinate, u = value of mean wall shear stress. The law of the wall constants are 0 = 0.41 and C+ = 5.7. Thus, Fig. 8 shows that the mean velocity profile obtained from PHASTA DNS is consistent with the law of the wall in the range of its applicability (y+ > 30).

1

10

y+

10 0

Fig. 8. PHASTA channel flow simulation velocity profile in wall coordinates (solid line) compared with the law of the wall (dashed line).

The computed values of interest have been used in the next step as the inflow boundary conditions for the NPHASE-CMFD RANS code. Note that in this test problem, the DNS-PHASTA case used periodic boundary conditions at the inflow and outflow of the domain. That means that the solutions at those planes are the same. Thus, taking the fluctuation measurements in the fully developed inflow results in the same evolutions as at the outflow. The left side of Fig. 6 shows how the fluctuating velocity field in the PHASTA solution has been averaged to provide a consistent mean velocity inflow for the NPHASE problem shown on the right side of Fig. 6. Fig. 9 shows the mean velocity and turbulent kinetic energy profiles at the inflow and outflow of the NPHASE domain. Note that the inflow represents the averaged PHASTA solution which is also shown in this figure. As can be seen, the High-Re k–ε model agrees well with the DNS results. Another test was also performed using the NPHASE-CMFD code by replacing the High-Re k–ε by a Low-Re version of the k–ε model. For this, more complex, problem a finer non-uniform mesh has been created, in which the first point off the wall was located at y+ = 1, and the overall grid has 50 cells across the channel. The averaging was performed over 30 periodic cycles (one cycle transfers a Lagrangian point through the domain in the stream-wise direction). Since it has been demonstrated by Yoder & Georgiadis (1999) that the Low-Re model results are nearly independent for wall mesh resolutions for y+ between 1, 2 and 5, the current wall resolution of y+ = 1 is appropriate. The averaged results obtained from PHASTA can be validated against a well-known Colebrook equation. We have obtained the liquid flow rate value of U = 1.008856 m/s in the PHASTA simulation by integrating the averaged velocity profile shown in Fig. 11. The periodic boundary conditions used in the PHASTA channel flow simulation require a pressure gradient as a driving force. The value of the pressure gradient used in the simulation is:

142

I.A. Bolotnov et al. / Nuclear Engineering and Design 247 (2012) 136–146

Fig. 9. Mean velocity, turbulent kinetic energy and turbulent dissipation rate profiles at the outflow of the PHASTA domain (averaged) compared to the outflow of NPHASE solution using high Reynolds number k–ε model.

Fig. 10. PHASTA turbulent channel flow domain (left) outflow data is time-averaged to obtain inflow boundary conditions for NPHASE (Low Reynolds number k–ε model case) domain (right).

∇ p = 0.004244 Pa/m. Thus, the friction factor observed in the presented DNS simulation can be estimated as: fDNS =

∇ pDh = 0.0333 U 2 /2

(20)

where Dh = 4.0 m and  = 1.0 kg/m3 . We can use Colebrook equation to estimate the friction factor analytically:

  1  = 2.035 log 0.25 Reh f + 0.142

(21)

f

where Reh = UDh / = 11, 033 is the Reynolds number based on the hydraulic diameter,  = 0.00036574 kg/ms is the liquid dynamic viscosity. The value of the friction factor computed from the Eq. (21) is equal to 0.03168 which compares well with the DNS observed value computed in Eq. (20). Fig. 10 shows how the NPHASE-CMFD mean velocity profile in the fine-grid case correlates with the PHASTA fluctuating velocity

field. Fig. 11 compares the mean velocity profiles of the averaged PHASTA-DNS result and NPHASE outflow solution as well as turbulent kinetic energy and turbulent dissipation rate distribution. Good agreement has been observed. The ability to process the DNS results to evaluate and extract the: mean velocity, turbulent kinetic energy and turbulent dissipation rate, for arbitrary geometries is crucial for the linking of the PHASTA-DNS code to the NPHASE-RANS code. The test cases discussed above provide an excellent testing ground for this purpose. The results shown in this section confirm that both k–ε models implemented in the NPHASE-CMFD code properly predict the parameters characterizing single-phase flow turbulence. Whereas the applicability of similar models to two-phase flows has been already confirmed for selected situations (Bolotnov et al., 2011), their use for specific situations and flow conditions should be approached with caution. Selected results obtained for the SFR accident scenario discussed in Section 2 are shown next.

Fig. 11. Mean velocity, turbulent kinetic energy and turbulent dissipation rate profiles of the time-averaged PHASTA outflow used as an inflow for the NPHASE problem compared with the NPHASE outflow result.

I.A. Bolotnov et al. / Nuclear Engineering and Design 247 (2012) 136–146

143

Fig. 12. A 2-D NPHASE simulation of side jet problem: (a) problem geometry, (b) stream-wise velocity distribution, and (c) gas volume fraction distribution.

6. NPHASE simulations of gas jet injection A two-dimensional analysis for a side jet problem was initially performed, which physically corresponds to the case of an elongated circumferential cladding breach. Naturally, turbulence was

still treated as a 3-D phenomenon. Fig. 12(a) shows the domain’s dimensions. As can be seen, the gas jet enters a fully developed turbulent flow of reactor coolant. Since in the present case, the gas is injected via a small opening (breach) in the fuel element cladding, it forms small sub-mm diameter bubbles, dispersed in the liquid

Fig. 13. 2-D NPHASE simulation of the side jet problem. Gas volume fraction (light green triangles), mean liquid velocity (blue diamonds), gas velocity (purple crosses) and turbulent kinetic energy (red squares) profiles are shown at different distances from the inflow. The left axes are for liquid/gas velocity and volume fraction plots, the right axes are for turbulent kinetic energy. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

144

I.A. Bolotnov et al. / Nuclear Engineering and Design 247 (2012) 136–146

Fig. 14. Coolant inflow velocity distribution.

metal coolant. In the presence of gravity, the bubbles accelerate in the flow direction, pulling the fluid around them thanks to the effect of drag. This results in a local increase of the liquid velocity, as shown in Fig. 12(b) and (c). As it is shown in Fig. 13, the local acceleration of liquid leads to higher velocity gradients which, in turn, results in the production of turbulence and an increase of the turbulent kinetic energy. Note that in the presence of bubbles, the relative velocity between the gas and liquid phases is nearly constant. A gradual spreading of the gas volume fraction profile when moving away from the injection point can also be observed in Fig. 13. In the final series of tests, the NPHASE-CMFD code has been used to simulate complete 3-D gas jet injection into the turbulent coolant

flow. The computational domain used in the computations is shown in Fig. 4. Details of the corresponding cross-sectional grid are shown on Fig. 14. As can be seen, the computational grid properly reflects the effect of geometry on the anticipated flow conditions, such as the use of multiple layers parallel to the fuel element walls to capture the impact of boundary layer on the flow field in this region. It should also be mentioned here that an issue which has been addressed in all studies involving the NPHASE-CMFD-based simulations deals with the impact of discretizations schemes on the results of calculations. In fact, assuring that the results are not affected by the mesh size or geometry is an important element of the standard Best Practice Guidelines. Details of addressing the mesh-independence issue, as well as several other aspects of the mathematical and physical consistency of model formulation in NPHASE-CMFD, and of the accuracy of predictions for a variety of experimental conditions, can be found in the references given in Section 4. Since NPHASE-CMFD uses the RANS approach to turbulence, which deals with averaged quantities rather than fluctuations, it can be assumed that the jet shape is symmetrical and thus, half of the domain can be used. It can be seen from the results presented on Figs. 15–18 that the domain size is adequately chosen to study the two-phase jet development in the 3-D domain. Parametric calculations have been performed to study the effects of different forces on the void fraction distribution. In the first series of runs, the effect was analyzed of bubble diameter in the dispersed gas phase. It was assumed in those runs that a side jet of fission gas (nitrogen properties) enters a reactor coolant channel filled with liquid sodium at turbulent flow conditions. A fully developed profile for the coolant inflow has been obtained using periodic boundary conditions corresponding to a single-phase flow. For the accident scenario discussed in Section 4, the mean jet velocity was 5.0 m/s, and the mean coolant velocity was 1.0 m/s. The mesh contained 18,000 elements, and the simulation was performed using a

Fig. 15. Volume fraction distribution for a jet with 0.125 mm bubbles diameter.

Fig. 16. Volume fraction distribution for a jet with 0.25 mm bubbles diameter.

I.A. Bolotnov et al. / Nuclear Engineering and Design 247 (2012) 136–146

145

Fig. 17. Volume fraction distribution for a jet with 0.5 mm bubbles diameter.

Fig. 18. Volume fraction distribution for a jet with 0.75 mm bubbles diameter.

High-Re k–ε model (Jones and Launder, 1972). The jet inflow velocity considered herein corresponds to a pressure of fission gas of 0.1 kPa. The jet velocity corresponds to the Mach number of 0.008 (considering the speed of sound of 625 m/s at 700◦ C). Thus, the compressibility effects did not have to be considered. Three different values of bubble diameter was used in the individual runs: 0.125 mm, 0.25 mm and 0.5 mm. The predicted gas concentration for the 0.125 mm bubbles is shown in Fig. 15. Note that this Figure has been compressed by a factor of 5 in the stream-wise direction to show details in the other directions. It can be noticed that a high bubble concentration near the injection point gradually decreases as the bubbles get dispersed through the domain and spread over an increasing fraction of the channel cross-section. Similar results for the 0.25 mm, 0.5 mm and 0.75 mm bubbles are shown in Figs. 16–18, respectively. As it can be seen, an increase in

bubble diameter results in a noticeable gas bubble departure from the walls. The stream-wise evolution of gas volume fraction for 0.5 mm bubbles, as well as the corresponding liquid and gas velocities are shown in Fig. 19. As expected, a rapid (nearly exponential) decrease in gas volume fraction along the flow is observed. We also can see that the effect of the buoyancy force on the relative velocity between the two phases. Fig. 20 provides a cross stream view of the distribution of void fraction near the domain outlet. A comparison between the cases with various bubble diameters reveals the fact that at practically any location inside the channel, the local void fraction decreases with the increasing bubble diameter. This is consistent with the fact that the drag force is inversely proportional to bubble diameter and, thus, the effect of phasic slip diminishes with decreasing bubble size so that the void fraction

Fig. 19. Stream wise evolution of gas volume fraction (dashed), coolant velocity (solid) and fission gas velocity (dotted) at the centerline location for a jet with 0.5 mm bubbles diameter. The horizontal axis indicated the distance from the coolant inflow.

Fig. 20. Cross stream distribution of void fraction (dashed), coolant velocity (solid) and gas velocity (dotted) near domain outlet for a jet with 0.5 mm bubbles diameter. The horizontal axis indicated the distance from the wall with the jet.

146

I.A. Bolotnov et al. / Nuclear Engineering and Design 247 (2012) 136–146

gradually approaches that limit corresponding to zero-interfacial slip conditions. 7. Conclusions It has been demonstrated that a DNS model can be used to obtain consistent boundary conditions for a RANS model of turbulence. The latter model was based on the modeling framework that has been extensively tested and validated before for both single and multiphase flows in systems of various geometries and operating conditions. The results presented in this paper, and the associated parametric testing of both physical and numerical assumptions, provide additional evidence that the overall NPHASE-CMFD-based computational model is numerically accurate and robust. The combined results of: (a) testing the physical models of turbulence and, (b) of the analysis on the effect of bubble size on the calculated velocity and phase distribution during fission gas injection into reactor coolant channels following local cladding failure in SFR, further indicate that the proposed model is capable of capturing three-dimensional phenomena governing local gas/liquid interactions during reactor accidents. The presented approach will be expanded in the future to include the phenomena governing cladding breach, detailed scenarios of multicomponent fission gas discharge, multiple rod geometry of the SFR fuel assembly, and others. Acknowledgments The authors would like to acknowledge the financial support provide to this study by the U.S. Department of Energy and the computation resources provided by the Computational Center for Nanotechnology Innovations (CCNI) and Scientific Computational Research Center (SCOREC) at Rensselaer Polytechnic Institute (RPI). References Anglart, H., Nylund, O., Kurul, N., Podowski, M.Z., 1997. CFD simulation of flow and phase distribution in fuel assemblies with spacers. Nuclear Science and Engineering 177, 215–228. Antal, S., Ettorre, S., Kunz, R., Podowski, M.Z., 2000. Development of a next generation computer code for the prediction of multicomponent multiphase flows. In: Proceeding of the International Meeting on Trends in Numerical and Physical Modeling for Industrial Multiphase Flow, Cargese, France. Antal, S.P., Podowski, M.Z., 2003. The effect of steam/water properties feedback on flow oscillations in boiling water reactors. In: Proceedings of the Tenth International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-10), Seoul, Korea. Antal, S.P., Podowski, M.Z., Lahey, R.T., Barber, D., Delfino, C., 2005. Multidimensional modeling of developing two-phase flows in a large adiabatic riser channel. In: Proceedings of the Eleventh International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-11), Avignon, France. Bertodano, M., Lee, S.J., Lahey Jr., R.T., Drew, D.A., 1990. The prediction of two-phase turbulence and phase distribution phenomena using a reynolds stress model. Journal of Fluids Engineering 112, 107–113. Bolotnov, I.A., Behafarid, F., Shaver, D., Antal, S.P., Jansen, K.E., Samulyak, R., Wei, H., Podowski, M.Z., 2010. Multiscale computer simulation of fission gas discharge during loss-of-flow accident in sodium fast reactor. In: Proceedings of the OECD Nuclear Energy Agency & IAEA Workshop (CFD4NRS-3) on Experimental Validation and Application of CFD and CMFD Codes to Nuclear Reactor Safety Issues, Washington, DC, USA. Bolotnov, I.A., Jansen, K.E., Drew, D.A., Oberai, A.A., Lahey Jr., R.T., Podowski, M.Z., 2011. Detached direct numerical simulation of turbulent two-phase bubbly channel flow. International Journal of Multiphase Flow 37, 647–659. Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surface tension. Journal of Computational Physics 100 (2), 335–354.

Chien, K.-Y., 1982. Predictions of channel and boundary-layer flows with a LowReynolds-number turbulence model. AIAA Journal 20 (1), 33–38. Drew, D., Passman, S.L., 1999. Theory of Multicomponent Fluids. Springer-Verlag, NY. Franca, L.P., Frey, S.L., 1992. Stabilized finite-element methods II; the incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics and Engineering 99 (2–3), 209–233. Gallaway, T., Antal, S.P., Podowski, M.Z., Guillen, D., 2007. NPHASE predictions of turbulent flow in the lower plenum of a prismatic gas-cooled reactor. In: Proceedings of the Twelfth International Topical Meeting on Nuclear Reactor Thermal-Hydraulics (NURETH-12), Pittsburgh, PA. Gallaway, T., Antal, S.P., Podowski, M.Z., 2008. Multidimensional model of fluid flow and heat transfer in Generation-IV supercritical water reactors. Nuclear Engineering and Design 238, 1909–1916. Jansen, K., 1999. A stabilized finite element method for computing turbulence. Computer Methods in Applied Mechanics and Engineering 174, 299–317. Jansen, K.E., 1993. Unstructured Grid Large Eddy Simulations of Wall Bounded Flows. Annual Research Briefs, Center for Turbulence Research, NASA Ames/Stanford University, 151. Jones, W., Launder, B., 1972. The prediction of laminarization with a two-equation model of turbulence. International Journal of Heat and Mass Transfer 15, 301–314. Laufer, J., 1954. The Structure of Turbulence in Fully Developed Pipe Flow. U.S. National Advisory Committee for Aeronautics (NACA), Technical Report 1174. Moser, R., Kim, J., Mansour, N., 1999. Direct numerical simulation of turbulent channel flow up to Re␶ =590. Physics of Fluids 11, 943–945. Nagrath, S., Antal, S.P., Podowski, M.Z., 2001. Multidimensional simulations of twophase flows in large volumes with injection spargers. In: Proceedings of the 4th International Conference on Multiphase Flow, ICMF’2001, New Orleans, LA, USA. Nagrath, S., Jansen, K.E., Lahey, R.T., Akhatov, I., 2006. Hydrodynamic simulation of air bubble implosion using a FEM based level set approach. Journal of Computational Physics 215, 98–132. Podowski, M.Z., 2009. On the consistency of mechanistic multidimensional modeling of gas/liquid two-phase flows. Nuclear Engineering and Design 239, 933–940. Pope, S., 2000. Turbulent Flows. Cambridge University Press. Sahni, O., Mueller, J., Jansen, K.E., Shephard, M.S., Taylor, C.A., 2006. Efficient anisotropic adaptive discretization of the cardiovascular system. Computer Methods in Applied Mechanics and Engineering 195 (41–43), 5634–5655. Sato, Y., Sekoguchi, K., 1975. Liquid velocity distribution in two-phase bubble flow. International Journal of Multiphase Flow 2, 79–95. Shaver, D.R., Antal, S.P., Podowski, M.Z., 2011. Multidimensional Modeling of Interfacial Boiling and Condensation. ANS Transactions, 2011 Winter Meeting, Washington, DC, USA. Taylor, C.A., Hughes, T.J.R., Zarins, C.K., 1998. Finite element modeling of blood flow in arteries. Computer Methods in Applied Mechanics and Engineering 158 (1–2), 155–196. Tejada-Martinez, A.E., Jansen, K.E., 2005. A parameter-free dynamic subgrid-scale model for large-eddy simulation. Computer Methods in Applied Mechanics and Engineering 194 (9), 1225–1248. Tiwari, P., Antal, S.P., Podowski, M.Z., 2006. Three-dimensional fluid mechanics of particulate two-phase flows in U-bend and helical conduits. Physics of Fluids 18 (4), 1–18. Tiwari, P., Antal, S.P., Podowski, M.Z., 2009. Modeling shear-induced diffusion force in particulate flows. Computers and Fluids 38, 727–737. Trofimova, A.V., Tejada-Martinez, A.E., Jansen, K.E., Lahey Jr., R.T., 2009. Direct numerical simulation of turbulent channel flows using a stabilized finite element method. Journal of Computers and Fluids 38 (4), 924–938. Tselischeva, E.A., Antal, S.P., Podowski, M.Z., Guillen, D., Beyer, M., Lucas, D., 2010a. Analysis of developing gas/liquid two-phase flows. In: Proceedings of the 7th International Conference on Multiphase Flow (ICMF-2010), Tampa, FL, USA. Tselischeva, E.A., Antal, S.P., Podowski, M.Z., 2010b. Mechanistic multidimensional analysis of horizontal two-phase flows. Nuclear Engineering and Design 240, 405–415. Whiting, C.H., Jansen, K.E., 2001. A stabilized finite element formulation for the incompressible Navier-Stokes equations using a hierarchical basis. International Journal of Numerical Methods in Fluids 35, 93–116. Wierzbicki, B., Antal, S.P., Podowski, M.Z., 2007. On the modeling of gas bubble evolution and transport using coupled level-set/CFD method. Nuclear Technology 158 (2), 261–274. Yoder, D.A., Georgiadis, N.J., 1999. Implementation and validation of the chien k–␧ turbulence model in the wind Navier-Stokes code. In: 37th Aerospace Sciences Meeting & Exhibit, Reno, Nevada, USA.