- Email: [email protected]
Simulation scaling studies of reactor core two-phase flow using direct numerical simulation
Nuclear Engineering and Design 358 (2020) 110435
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
Simulation scaling studies of reactor core two-phase flow using direct numerical simulation
T
⁎
Joseph J. Cambareria, , Jun Fangb, Igor A. Bolotnova a b
Department of Nuclear Engineering, North Carolina State University, Raleigh, NC 27695, United States Nuclear Science and Engineering Division, Argonne National Laboratory, Lemont, IL 60439, United States
A R T I C LE I N FO
A B S T R A C T
Keywords: DNS Interface tracking PWR subchannel Scaling studies
Tremendous growth in supercomputing power in recent years has resulted in the emergence of high-resolution flow analysis methods as an advanced research tool to evaluate single and two-phase flow behavior. In particular, unstructured mesh-based methods have been applied to analyze flows in complex reactor core geometries, including those of light water reactors (LWR). The finite-element based code, PHASTA, is utilized to perform large-scale simulations of two-phase bubbly flows in LWR geometries. Given the large computational cost of direct numerical simulation (DNS) coupled with interface tracking methods (ITM), typical domains encompass a portion of a single subchannel. In the presented research, the state-of-the-art analysis of turbulent two-phase flows in complex LWR subchannel geometries are demonstrated at both prototypical reactor parameters as well as scaled low pressure conditions. Three different cases are studied, a high-pressure simulation in prototypical reactor subchannel geometry, a low-pressure case in prototypical geometry and a final low-pressure case in a geometry scaled up to conserve the ratio between the bubble size and the domain pitch. Utilizing advanced statistical processing tools, these simulation conditions are compared to shed light on the relevancy of two-phase flow characteristics given the significant differences between LWR and low-pressure conditions. These findings can lead to the generation of useful guiding principles when researchers need to scale the two-phase flow behavior captured at low pressure and temperature conditions to those at reactor operating conditions.
1. Introduction Due to the high heat transfer coefficients observed within the subcooled boiling regime, the latter is one of the most efficient forms of convective heat transfer. This is in part due to the turbulence intensity generated by the cyclic generation and condensation of vapor bubbles within the flow. However, this regime is limited by the critical heat flux of the system, where past this point highly reduced heat transfer coefficients are observed, resulting in surface burnout or fuel rod damage in nuclear reactor cores. Consideration of the critical heat flux is particularly important in the thermal–hydraulic (TH) and safety analysis of nuclear reactor systems to maintain fuel rod integrity during normal operation as well as accident scenarios. With this concern in mind, the accurate prediction and modeling of turbulent two-phase flows is desired to guide LWR design (Turinsky and Kothe, 2016). This drives the need for expanding the understanding of the physical processes and mechanisms behind turbulent two-phase flows at both high- and lowpressure conditions. For this purpose, multiple TH codes have been developed to support
⁎
the efficient and safe design of nuclear reactors, such as RELAP (Fletcher and Schultz, 1995), TRACE (Odar et al., 2003) and COBRA-TF (Gluck, 2008). While highly cost efficient, the drawback of these TH codes is that they still rely on 1-D or 0-D two-phase closure models. This reliance means that the predictive capabilities of these thermal-hydraulics codes are heavily dependent on the accuracy and applicability of the applied closure relations. Due to the extreme and complex nature of LWR operating parameters, it is highly challenging and expensive to perform high-resolution (interface capturing) experimental measurements at pressures and temperatures observed in nuclear reactors. This results in a lack of quality experimental data for high fidelity model development. For this reason, many subcooled boiling models rely on experimental data from simplified geometries, such as that of Kroeger and Zuber (1968) or the scaling of low-pressure refrigerants to the high pressure and temperature conditions typical of reactor operation. However, this introduces new complexities such as the applicability of scaling relations between low-pressure and high-pressure turbulent bubbly flows and the restrictions of the geometry studied. The breakdown and lack of scalability of these models to low pressure subcooled
Corresponding author. E-mail addresses: [email protected] (J.J. Cambareri), [email protected] (J. Fang), [email protected] (I.A. Bolotnov).
https://doi.org/10.1016/j.nucengdes.2019.110435 Received 30 May 2019; Received in revised form 22 October 2019; Accepted 11 November 2019 0029-5493/ © 2019 Elsevier B.V. All rights reserved.
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
adequately handle turbulence and bubble interactions with and around complex geometries, such as spacer grids or mixing vanes within a LWR subchannel (Fang et al., 2018). Introducing the novel bubble tracking method, which was recently developed, allows for the collection of numerical data on the individual bubble level for the use in advanced statistical post processing (Fang et al., 2013). With the ability to have highly scalable performance on massively parallel systems, up to 786x1024 processors on the IBM Blue Gene/Q Mira system at Argonne National Laboratory, PHASTA is a viable research tool for the study and advanced modeling of turbulent two-phase flows (Rasquin et al., 2014). Mesh generation employs the use of Chef, a derivative of the open source parallel unstructured meshing infrastructure (PUMI) (Ibanez et al., 2015), as well as utilizing tools from Simmetrix, Inc. PHASTA, coupled with Chef meshing tools, is able to work with hexahedral, tetrahedral and mixed finite element meshes along with being able to perform adaptive mesh refinement and coarsening.
boiling conditions for water have also been reported, where unsatisfactory predictions of void fraction distribution are observed. Woodruff et al. (1997) observed this behavior when testing RELAP5 against experimental transient data for the SPERTIV-D-12/25 test. There have also been efforts to improve the modeling of subcooled boiling flows. Hari and Hassan (2002) worked to enhance the predictive capabilities of RELAP5 by implementing the heat partitioning model of Zeitoun and Shoukri (1997) and modifying the interfacial condensation rate. The deficiencies of subcooled boiling models for computational fluids dynamics (CFD) have also been reported by Tu and Yeoh (2002) and Končar et al. (2004), who each proposed multidimensional twofluid subcooled boiling models for the general purpose CFD code, ANSYS-CFX. This indicates a need for further investigation to create guiding principles for researchers when scaling between low-pressure conditions and conditions typical to reactor operation. Generally, better closure relations are required to define the wall heat flux partitioning, the bubble departure diameter and bubble dynamics within the flow. Recently, the expansion of high performance computing (HPC) facilities has led to the viability of utilizing DNS to supplement experimental data when developing accurate two-phase closure models for advanced reactor design (Fang et al., 2018; Fang and Bolotnov, 2017). While DNS is still one of the most computationally expensive methods, the achievable scale and mesh resolution of DNS has observed tremendous growth. Given adequate spatial and temporal resolution, DNS can represent all relevant scales of turbulence without the incorporation of any turbulence models. Coupled with interface tracking methods, DNS can provide high resolution numerical data for turbulent single- and two-phase flows (Moser et al., 1999; Bolotnov, 2013). This ability to collect reliable numerical data has led to DNS becoming an accepted source to supplement experimental data. This work focuses on the study of how the system pressure affects bubble dynamics and interactions with core structures, such as spacer grids and mixing vanes, within LWR geometries. The LWR subchannel geometry is selected to analyze these effects due to the importance of understanding bubble behavior near the mixing vane region and its importance in the calculation of heat transfer mechanisms during normal reactor operation and accident scenarios. Three different simulations are considered. The first recreates the conditions that are expected for typical reactor operation within an LWR subchannel geometry and the analysis of this case has been previously presented in literature (Fang et al., 2018). The second and third simulations scale the case parameters to those expected at atmospheric conditions, with one simulation conserving the geometry and the other scaling the geometry to match the new bubble diameter to hydraulic diameter ratio. In both simulations, the bubble deformability, through the Eötvös number (Eo), are preserved by increasing the bubble diameter with the differing fluid properties. For each simulation, a previously developed bubble tracking method was applied to track the statistics of individual bubbles within the flow (Fang et al., 2013). These new simulations and state-of-the-art statistical analysis of turbulent two-phase flows at both low-pressure conditions and prototypical reactor operating conditions will help to shed light onto relevant bubbly flow behavior in LWR geometries and the impact of spacer grid and mixing vanes on two-phase models and closure laws.
2.1. Flow solver The conditions considered within this research describe an isothermal and incompressible system. PHASTA solves the weak form of the Incompressible Navier-Stokes (INS) equations in three dimensions through a stabilized FEM. The spatial and temporal discretization for the INS equations within PHASTA have been previously discussed at length by Whiting and Jansen (2001) as well as Nagrath et al. (2005). The governing equations for continuity and momentum in the strong form of the INS equations are given by:
∂ui =0 ∂x i
ρ
∂τij ∂p ∂ui ∂u + fi + + ρuj i = − ∂x j ∂x i ∂t ∂x j
(1)
(2)
here, ui denotes the velocity for the i-th dimension with i = 1, 2 and 3, ρ the fluid density, p the static pressure, τij as the viscous stress tensor and finally fi represents the body force in the i-th coordinate, including the gravitational and surface tension forces. This final term includes both the gravitational and surface tension forces. The viscous stress tensor for an incompressible flow of a Newtonian fluid is connected to the dynamic viscosity, μ , of the fluid and the strain rate tensor, Sij .
∂uj ⎞ ∂u τij = 2μSij = μ ⎜⎛ i + ⎟ ∂x i ⎠ ∂ x j ⎝
(3)
Expanding from single- to two-phase capabilities, PHASTA employs the level-set method to distinguish between the liquid and gas phase (Sussman et al., 1998, 1999). The level-set method creates a precise, continuous, signed distance function, φ , which represents the bubble interface at the zero of the level-set distance field. Thus, φ = 0 represents the interfacial boundary between the phases. By convention, a positive value of the level-set distance field indicates the liquid phase, while the gas phase is denoted by a negative value. The scalar, φ , is advected with the fluid according to the advection equation.
∂φ ∂φ =0 + ui ∂x i ∂t
2. Numerical methods
(4)
To account for the discontinuity created by the fluid properties at the interface, a smoothed Heaviside kernel function, H I µ̂ , is introduced (Sussman et al., 1999). This function assumes a finite interfacial thickness, and the properties over the interface, the density and viscosity of each phase, are smoothed to remove the jump discontinuity, as presented below.
The finite element method (FEM) research code, PHASTA (Jansen, 1999), was used to perform the direct numerical simulations. PHASTA is a parallel, hierarchic, higher-order accurate, adaptive, stabilized transient analysis flow solver for the solution to both incompressible and compressible flows. The flow solver has been shown to be a highly effective research tool for a broad range of turbulent flow conditions and scales, having been applied to RANS, LES, Detached Eddy Simulation (DES) and DNS cases (Jansen et al., 1993; Whiting and Jansen, 2001). With its ability to process unstructured grids, PHASTA can 2
ρ (φ) = ρ1 Hε (φ) + ρ2 (1 − Hε (φ))
(5)
μ (φ) = μ1 Hε (φ) + μ 2 (1 − Hε (φ))
(6)
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
2.2. Bubble tracking method Throughout the history of DNS, there have been multiple techniques for numerical data collection and analysis. Single- and two-phase quantities had been collected through the implementation of static virtual probes arrayed in 2-D planes within the computational domain. The virtual probes planes are placed normal to the streamwise direction, with probe locations set at specific distances from the wall. For each distance to the wall value there is a specific number of virtual probes where similar statistical behavior is expected to occur. Over the course of the simulation, these probes collect instantaneous numerical data, such as the phase velocities or void fraction. Statistical averaging in time and space is then used to compute two-phase turbulent flow parameters, such as the mean flow velocity, void fraction distributions, Reynolds stresses, turbulent kinetic energy or the turbulence dissipation rate of the liquid phase. While the virtual probe statistical analysis approach has made major contributions to the development and study of turbulent closure models (Bolotnov, 2013), especially for singlephase turbulent flows, the efficient and useful numerical data collection for two-phase turbulent flows presents more of a challenge. The virtual probe arrays are in static planar locations, removing their ability to adequately capture the spatial movement patterns and behavior of bubbles in two-phase flows, especially in complex geometries such as those with spacer grids and mixing vanes. For this reason, the development and implementation of a bubble tracking method was sought after. Developed by Fang et al. (2013), a scalar marker field is initialized along with the level set distance field. This scalar field individually tracks and collects the instantaneous numerical data of each separate bubble within the computational domain. Each bubble is assigned a unique ID, which is advected with the level set to dynamically track the bubbles simulated. Additionally, numerical data from the surrounding liquid shell is collected to assess the local flow conditions, such as the liquid shear or relative velocity, which the bubble is encountering. The bubble information that can be extracted from this method is highly valuable for the generation and improvement of models and closure laws, such as the drag or lift coefficient, for use in lower resolution computational fluid dynamics (CFD) to better predict turbulent bubble flows.
a)
b) Fig. 1. a) Computational domain setup for the three cases (front walls hidden to display inner structures) and b) location of the measurement planes within the domain.
spacer grid and mixing vanes, as displayed in Fig. 1 below. Here, the front walls of the domain are hidden to display the inner structures of the spacer grid and mixing vanes. The length of the subchannel domain is 4.43 rod diameters in length (40.5 mm). This domain length was selected to balance the quality of numerical statistics and the computational cost of the simulations. Two measurement planes of static virtual probes are located within the computational domain, set at locations of 0.86 and 2.28 rod diameters from the end of the mixing vanes (2.74 and 4.16 rod diameters from the domain origin), as displayed by the red dashed lines in Fig. 1. Three different simulation setups are considered in this investigation: a 1% void fraction (262 bubbles) highpressure case previously presented in literature (Fang et al., 2018, 2019), a simulation that conserves the prototypical subchannel geometry while scaling to low-pressure conditions and a final simulation which scales the geometry size to conserve the bubble diameter to pitch ratio for low-pressure conditions. As a note, for the following comparison the simulations will be referred to as HP for the high-pressure simulation, LP1 for the case conserving the initial geometry and scaling to low-pressure conditions and LP2 for the simulation which scales both the pressure and geometric conditions. Periodic boundary conditions are applied at the inlet and outlet of the domain as well as the transverse faces to allow for the turbulence and bubble interactions with adjacent subchannels without the need for meshing those regions. Along the surface of the fuel rods, spacer grid and mixing vanes the no-slip wall boundary condition is applied.
3. Simulation setup When considering the thermal–hydraulic analysis of an LWR, the subchannel geometry plays a critical role in the overall reactor core design. As coolant flows through the subchannel, it will carry heat generated from the reactor fuel rods to the steam generators. During prototypical PWR operating conditions, subcooled boiling can occur in the upper region of the reactor core, resulting in the production of steam bubbles within the subchannel coolant. The inclusion of vapor bubbles within the flow typically results in an increase of heat transfer efficiency because of the bubble induced turbulence and the subcooled boiling phenomenon. On top of this, the inclusion of complex geometries within the subchannel, such as spacer grids and mixing vanes, can greatly affect the solution of the coolant flow field, enhancing the turbulent intensity of the coolant and increasing the heat transfer and heat removal from the fuel rod surfaces. While both vapor production and complex geometries further complicate the heat transfer mechanisms at play, the full impact these effects on two-phase flows within LWR subchannels are not yet fully understood due to the difficulty of obtaining high-resolution experimental measurements at the extreme temperatures and pressures observed in the reactor core environment.
3.2. Parameter scaling to low pressure simulations Detailed results of the 1% void fraction case at conditions approaching prototypical reactor operating conditions have been presented in literature (Fang et al., 2018, 2019). For this simulation, the fluid properties were set by using the saturated properties of water and vapor at 573 K. The saturation temperature selected for this work is associate with a corresponding system pressure of 8.58 MPa, displayed in Table 1, which is equivalent to conditions observed within an LWR. However, in the push to create guiding principles for researchers to follow when scaling between pressure conditions, the following research focuses on the comparison between the numerical data gathered from those investigations to DNS of turbulent bubble flows completed at low-pressure conditions. For the low-pressure simulations, the system pressure was selected to be 0.1 MPa and the fluid properties are based
3.1. Computational domain design Given the available computational resources, and knowing the importance of the subchannel geometry, the computational domain generated for this work includes a single subchannel with a reduced size 3
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
Table 1 Property comparison at prototypical reactor conditions and low-pressure conditions. Parameter/Case
High-pressure
Low-pressure
Unit
Pressure Saturation Temperature Fluid Density (liquid; vapor) Fluid Dynamic Viscosity (liquid; vapor) Surface Tension Coefficient Gravitational Acceleration Bubble Diameter Eötvös Number
8.58 573 712; 46 8.59E−05; 1.97E−05 0.014 9.81 0.65 0.25
0.10 373 959; 0.59 2.82E−04; 1.23E−05 0.059 9.81 1.25 0.25
MPa K kg/m3 Pa s N/m m/s2 mm –
Table 2 Overview of significant geometric and computational properties. Parameter/Case
HP
LPC1
LPC2
Bubble diameter (mm) Hinze Scale (mm) Rod diameter (mm) Rod pitch (mm) Domain height (mm) Mesh cell size in bulk region (μm) Mesh cell size of the first boundary layer (μm) Total number of mesh cells Number of bubbles resolved Mean velocity (m/s) Flow Reynolds number
0.65 0.82 9.14 12.60 40.50 29.30 8.14
1.25 1.57 9.14 12.60 40.50 57.21 7.15
1.25 2.04 17.66 24.33 78.22 61.68 7.71
1.10 billion 262 0.75 8.1E04
0.29 billion 37 0.75 3.3E04
1.13 billion 262 0.75 6.4E04
Fig. 3. Scaled virtual probe locations for the LP2 simulation.
of the bubbles had to be approximately doubled to ensure the same value of the Eötvös number. By guaranteeing the conservation of this dimensionless characteristic, the behavior of the bubbles in the prototypical case will be most readily compared to those in the low-pressure simulations.
on the corresponding saturated properties of water and vapor. Of main interest in this study was the comparison of bubble behavior and dynamics around complex LWR geometries (e.g. spacer grids and mixing vanes) when scaling to low pressure conditions. For this reason, the bubble deformability, through the Eötvös number, presented in (7) below, was conserved. Here, the Eötvös number is calculated using the equivalent diameter, db , of a perfectly spherical bubble. From the previous studies of turbulent bubbly flows in LWR subchannels, the resulting average Eötvös number of the resolved bubbles was 0.25. Because of the large change in simulation properties resulting from the change in pressure, mainly due to the influence of the larger surface tension coefficient at atmospheric pressures, the diameter
a)
Eo =
Δρgdb2 σ
(7)
It should be noted that the Reynolds number of the original case was not conserved throughout each simulation. At higher Reynolds numbers the bubbles are more likely to breakup as the influence of the external forces upon the bubble overcome the surface tension forces. The impact of this can be generally classified with the calculation of the Hinze scale, given in (8) below (Garrett et al., 2000). Here, c is a model constant with a value of 0.363 and ε is the dissipation rate. If the bubble radius is greater than the Hinze scale it is more likely to breakup under the flow conditions. For the high- and low-pressure simulations the value of the Hinze scale are listed in Table 2 below. In each of these
b)
Fig. 2. Initial bubble size and geometry comparison of the high- and low-pressure cases for a) the conserved geometry and b) the scaled geometry. 4
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b) Fig. 4. Snapshots for the LP1, HP and LP2 simulations (left to right) at a) the two-phase initialization (0 BDT) and b) 30 BDT with the bubbles delineated by their IDs and the slices displaying the instantaneous velocity.
the turbulence of the flow. Since the both the high- and low-pressure bubble radius falls below the Hinze scale, the lack of Reynolds number scaling is not anticipated to have a large impact on the bubble
simulations the bubble size is much larger than the Kolmogorov length scale (on the order of 10 µm for these simulations) and smaller than the Hinze scale, so minimal breakup of the bubbles is to be expected due to 5
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
Fig. 5. The instantaneous liquid velocity field at the first measurement plane with the location of the outer tip of the mixing vanes marked (dashed black line).
a)
b)
Fig. 6. The liquid a) mean streamwise velocity and b) velocity tangential to the fuel rod surfaces at the first measurement plane location.
a)
b)
Fig. 7. The turbulent kinetic energy profile of the liquid phase versus the mixing vane tip location (black dashed line) for the a) first measurement plane and b) second measurement plane.
bubbles to only 37. To create a more direct comparison of the highpressure case, another low-pressure simulation was initialized, but the geometric properties of the high-pressure case were scaled to conserve the bubble diameter to hydraulic diameter ratio of the original geometry. To conserve this ratio, the new domain was scaled linearly by a factor of 1.93 times the size of the original domain. Additionally, to replicate the original high-pressure simulation conditions as close as possible, the initial bubble locations of the 1% void fraction case were scaled with the same geometric scaling factor. A comparison of
dynamics.
aH = c (σ / ρ)3/5ε −2/5
(8)
With the large change in the bubble diameter resulting from the conservation of the Eötvös number, the ratio of the bubble diameter to the geometric conditions of the original computational domain were greatly affected. To maintain the 1% void fraction initialization that was applied to the high-pressure simulation, the number of bubbles in the atmospheric pressure case had to be greatly reduced from 262 6
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b)
c) Fig. 8. Anisotropy distribution comparison for the a) HP and LP1 case, b) HP and LP2 case and c) the difference between the data sets where a one-to-one fit (dashed black line) is presented with lines for ± 20% (dotted black lines).
Now, the smaller turbulent eddies can be developed on the finer mesh for a more time efficient process. The refinement process can continue until the required turbulent scales are achieved. For each simulation, the desired turbulent profile and flow rate is accomplished by adjusting the flow driving pressure gradient in the streamwise direction until it balances the wall drag on the flow to achieve the desired flow rate. The inclusion of the spacer grid and mixing vanes creates a highly efficient process for the development of single-phase turbulence. The mean flow velocity is set at 0.75 m/s for each void fraction case, which corresponds to a hydraulic Reynolds number of 80,737 for the HP case, 33,097 for LP1 and 63,812 for LP2. The gravitational force is applied in the direction opposite to the flow to create the conditions of vertical upward channel flow. Once the single-phase turbulence reaches a statistically steady state, the two-phase conditions are initialized. Fig. 4 exhibits the initial bubble locations, colored by their ID numbers, for the HP simulation and the low-pressure cases, as well as displaying the initial instantons velocity field. Fig. 4 also displays a further two-phase time step for the simulations at approximately 30 bubble diameter time (BDT) based on the HP case. The BDT is a nondimensional unit that allows for a more direct comparison between the cases and is defined as the average time interval that a bubble required to move one diameter in distance. For the presented analysis the BDT of the LP1 and LP2 cases was converted to use the bubble size of the HP case. To ensure the quality of the statistics collected, the two-phase numerical data was collected once a quasi-steady state was reached for each simulation. High quality statistics are a necessity for numerical simulations and one of the easiest methods to accomplish this is to
important geometric and computational properties for the three cases are displayed in Table 2 below. To better visualize the change in the geometric and simulation parameters, Fig. 2 shows the high and low pressure bubble sizes with both the conserved and scaled geometries.
4. Results and discussion All of the cases were run on the IBM Blue Gene/Q Mira supercomputer at Argonne National Laboratory (ANL). The new LP1 and LP2 cases were carried out on a significantly larger computational allocation than previously awarded, so the mesh quality, based on the number of elements across a bubble diameter, was increased to better resolve the interfacial behavior. From the previous 1% void fraction case, the number of elements across a bubble diameter was increased from 21 to approximately 24. Both the conventional static array of probes and the novel bubble tracking algorithm data collection methods were applied for numerical data collection in the simulation scaling studies. Two measurement planes of static virtual probes are placed within the computational domain, set at locations of 0.86 and 2.28 rod diameters from the end of the mixing vanes for each case. For the LP2 case, the locations of the probes were scaled with the same geometric ratio as before. The resulting probe plane is presented in Fig. 3. The numerical data collected by the probe arrays is supplemented with that from the bubble tracking algorithm. For an efficient large-scale PHASTA simulation workflow, the large turbulent eddies are initially developed on a coarse mesh before the mesh is uniformly refined and the solution migrated onto a finer mesh. 7
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b)
c) Fig. 9. Anisotropy distribution of the off-diagonal components for the a) HP and LP1 case, b) HP and LP2 case and c) the difference between the data sets where a one-to-one fit (dashed black line) is presented with lines for ± 20% (dotted black lines). Table 3 Comparison of group bounds for the high- and low-pressure simulations. HP
LP1
LP2
Group Number
Spatial range (by wall distance, unit: mm)
Group Number
Spatial range (by wall distance, unit: mm)
Group Number
Spatial range (by wall distance, unit: mm)
F1 F2 F3 F4 F5 F6
(0.0, 0.51) (0.51, 0.98) (0.98, 1.41) (1.41, 1.82) (1.82, 2.47) (2.47, 4.34)
FR1
(0.0, 0.98)
FR2
(0.98, 1.82)
FR3
(1.82, 4.34)
F1 F2 F3 F4 F5 F6
(0.0, 0.99) (0.99, 1.89) (1.89, 2.71) (2.71, 3.51) (3.51, 4.77) (4.77, 8.69)
increase the sample size of the numerical data set. Due to the large computational cost and the transient nature of turbulent DNS simulations, this is not always an option though. To balance the quality of statistics gathered and the computational cost, special attention is given to the convergent behavior of the time averaged statistics. Three metrics were used to assess this condition; the convergence of the mean streamwise velocity, the balance of the bubble population in relation to the distance to the wall and the development of the bubble relative velocity. The numerical data collected from the virtual probe planes is averaged across multiple time windows, and the convergence of the streamwise velocity profile is evaluated. The second two methods employ the bubble tracking numerical data. The bubble population distribution and relative velocity can be dynamically tracked over the course of the entire simulation. The solution is assumed to reach a quasi-steady state condition when the relative vapor/liquid velocity has stabilized in the streamwise direction and the distribution of bubbles in
Fig. 10. The bubble tracking groups bounds for HP and LP2 (solid lines) and LP1 (dotted bands) in a quarter subchannel cross-section.
relation to the distance to the wall has stabilized. For this work, it was identified that approximately 40 BDT was required to reach this quasisteady state, at which point the streamwise velocity profile had also converged. Due to the fuel rod geometry, the local wall normal coordinates at the virtual probe locations change in reference to the global coordinate system. This presents a unique challenge in processing the numerical data, where the wall distance at two opposite fuel rod surfaces has a 8
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b)
Fig. 11. Evolution of the bubble population for the a) near wall and b) channel center distance to the fuel rod surface groups as a function of the distance from the domain inlet for the HP (solid lines), LP1 (dotted lines) and LP2 (dashed lines) cases.
from subchannels without the spacer grid and mixing vane do not show this trend, where a flatter profile is seen outside of the boundary layer region (Fang and Bolotnov, 2017). By the new coordinate system convention, the tangential direction is positive in a counterclockwise direction along the surface of the fuel rods. Due to the dominant clockwise rotation of the coolant, the tangential velocity shows a large negative velocity between the locations of 0.4 and 0.8 r/R. Within this region the mixing vanes have the largest influence upon the fluid flow. On top of the mean velocity profile, the static virtual probe planes also can process the velocity fluctuations in the domain to calculate the turbulent kinetic energy (TKE) profile for the liquid phase. The impact of the internal structures is also prominent on the turbulent kinetic energy profiles. Compared to a simple subchannel without an internal structures (Fang and Bolotnov, 2017), the presence of the spacer grid and mixing vanes significantly increases the magnitude of the TKE profile in the coolant flow, as seen in Fig. 7. Within the center of the subchannel, as well as in the near wall region for LP1, the low-pressure cases observe higher values for the TKE. However, for both of the low-
zero value. To account for this aspect of the subchannel geometry, the Cartesian coordinates within each quarter of the subchannel are adjusted into a streamwise, tangential and wall normal component system. To produce the time averaged statistics, the instantaneous numerical data from the virtual probes located at each distance to the wall value are averaged over a specific time window. For this study, the numerical data collected by the virtual probes within the rod gaps was excluded to reduce the influence of the cross flow on the statistics. With the inclusion of internal structures in the subchannel, especially resulting from the mixing vane geometry, the instantaneous liquid velocity field reveals a centrifugal pattern as observed in Fig. 5. The mixing vanes produce the rotating flow in the core of the subchannel, forcefully pushing the flow in the direction of the outer regions. This has the additional effect of generating a lateral flow at the surface of the fuel rods, enhancing cross flow and mixing between subchannels. The mean streamwise velocity, depicted in Fig. 6, also displays this effect, where a reduction in the liquid velocity in the center of the subchannel is observed. Previous mean velocity profiles 9
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b)
Fig. 12. The bubble relative velocity for the a) near wall and b) channel center distance to the fuel rod surface groups as a function of the distance from the domain inlet for the HP (solid lines), LP1 (dotted lines) and LP2 (dashed lines) cases.
prototypical reactor and low-pressure conditions using state-of-the-art statistical analysis based on the dynamic bubble tracking method. Through the implementation of the bubble tracking algorithm, detailed bubble information is collected over the entire lifetime of a bubble during the two-phase flow simulations. The collected bubble tracking data is divided into several bubble groups based on their distance to the fuel rod surface and the distance to the domain inlet for a more detailed study of the bubble dynamics involved. For this study of the HP case, the bubble numerical data was split into six groups based on their distance to the wall and fourteen groups based on the distance to the domain inlet. The group outer bounds were set such that each group has an equivalent cross sectional area, and the group bounds are presented in Table 3. Adjustments were also made for the group bounds in the LP1 case due to the size increase of the bubbles within the conserved subchannel geometry. Specifically, the number of distance to the fuel rod groups was reduced to three due to the size of the bubbles in relation to the geometry. A visual representation of the group bounds within a quarter of the subchannel cross section is displayed in Fig. 10 below. Within this figure, the solid lines with the dotted bands represent where the group bounds for the HP and LP2 cases coincide with that of the LP1 case. Since the initial bubble size is larger than that of the viscous sublayer for all cases, limited numerical bubble data can be gathered for the near wall relations. However, important bubble behavior and dynamics can still be evaluated within the turbulent log-law region of the subchannel geometry. To track the specific bubbles that occupy each group, the centroid of the bubble is recorded, and its location is compared in relation to the group bounds. While this is a fairly simple approach, future advancements of this technique can account for when
pressure cases the influence of the complex geometry fades away faster than that of the high-pressure case. Future studies will investigate the impact of bubble induced turbulence to isolate the dominate effect changing the TKE profile. With the application of DNS in this work, each term within the Reynolds stress tensor can be computed for the liquid phase. Specifically, the directional components of the Reynolds stress tensor −
are calculated for the streamwise (u'u' ), normal to the fuel rod surface − ' '
−
(v v ) and tangential to the fuel rod surface (w 'w ' ) components. The twophase flow anisotropy distribution can be computed by normalizing these values with the turbulent kinetic energy. For a flat plat boundary layer the ratio of the normal stresses is typically found to be −
−
−
u'u': v 'v ': w 'w ' = 8: 4: 6, so the subchannel anisotropic values are scaled by a factor of 9 for an easier comparison (Wilcox, 2002), however with the inclusion of the complex geometry in the subchannel domain this value is expected to differ. Fig. 8 presents the comparison of the highand low-pressure data as well as the difference between the values. The comparison in Fig. 8 presents all of the directional components of the Reynolds stress tensor in relation to a one-to-one and a difference of ± 20%. Overall, the LP2 case appears to more closely match the turbulent anisotropy distribution of the HP case. The off-diagonal components of the Reynolds stress tensor also play an important role in the overall characterization of the subchannel turbulence. Fig. 9 presents the comparison of the off-diagonal Reynolds stress tensor terms for the HP, LP1 and LP2 cases. While a wealth of information can be calculated from the turbulent statistics of DNS, this work mainly focuses on the comparison of 10
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b)
Fig. 13. The bubble deformability for the a) near wall and b) channel center distance to the fuel rod surface groups as a function of the distance from the domain inlet for the HP (solid lines), LP1 (dotted lines) and LP2 (dashed lines) cases.
Further analysis that can be completed with the bubble tracking method involve the tracking of novel statistics such as the bubble deformability factor. The bubble deformability factor, γd , is defined as the ratio of the minimum value of the level set field found within a bubble to the radius of a sphere with the equivalent volume as the considered bubble. In short, this term quantifies how deformed a bubble is, with an ideal sphere having a value of γd = 1.0 . Fig. 13 displays the effect that the complex geometry have upon the shape of bubbles within the flow, where bubbles entering and leaving the spacer grid and mixing vane region, especially in the core of the flow and close to the wall, are notably more deformed than those outside of the complex geometry region of influence. Within the spacer grid region this deformation is likely caused by flow confinement, and after the mixing vanes the bubbles tend to deform under the influence of the turbulent flow. Both of these effects decay as the bubbles move downstream from the region of influence of the complex geometry. For both low-pressure cases, the bubbles are significantly more deformed than that of the high-pressure case, with LP1 observing the highest amounts of deformation. This deformability is most likely a result of the relatively large bubble size moving though the complex geometry coupled with the increased weight of the low-pressure fluid. Additionally, for the central bubble groups, the bubbles tend to return to a more spherical shape after leaving the area of influence generated by the spacer grid and mixing vanes within the channel.
a bubble is moving across a group bound and indicate the bubbles that are entering and leaving a group. Separating the bubbles into groups based on the distance to the fuel rod allows for the evaluation of the population distribution and bubble relative velocity throughout the subchannel over the course of the bubbly flow simulations, as seen in Fig. 11 and Fig. 12. For this work the relative velocity convention is taken as the bubble velocity minus the liquid velocity. Within these plots the solid vertical black lines indicate the beginning and the end of the spacer grid region, while the dashed vertical black line specifies the end of the mixing vanes. This method allows the observation of this migration trend throughout the length of the subchannel and the evolution of bubble dynamics as they are affected by the spacer grid and mixing vanes. Generally, upward two-phase flows typically have a void fraction profile that peaks closer to the wall, and previous subchannel simulations excluding the spacer grid and mixing vanes demonstrate this migration trend of bubbles towards the wall (Fang and Bolotnov, 2017). However, the centrifugal effects of the mixing vanes alter the dynamics of the flow, forcing bubbles to migrate into the center of the subchannel after they encounter the mixing vanes. Notably, after the mixing vane region, the most central groups (F6 and FR3) see large increase in bubble population due to the turbulent mixing effects of the mixing vanes. This increase in population is accompanied by a sharp decrease in the relative velocity of the bubbles found within those groups, although this decrease is observed for every distance to the fuel rod surface group. Both the LP1 and LP2 cases observe a larger migration of bubbles into the core of the flow, however a smaller relative velocity drop when compared to the HP case.
5. Conclusion In the presented research, turbulent bubbly flows within LWR 11
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
SCOREC at RPI (Chef) and Kitware (ParaView).
subchannels are investigated at conditions prototypical to reactor operation and those found in low-pressure scenarios. Using high-fidelity, state-of-the-art interface tracking simulations, all relevant scales of liquid turbulence are resolved, and the interface is captured with the level-set method. Three different cases are considered to carry out the scaling study of the effects of system pressure upon bubble dynamics within the flow. Using static virtual planes, the impact of the complex geometry can be observed on both the mean flow velocities and the TKE profile. Further investigations into the effect of bubble induced turbulence and the change in liquid properties are planned to examine its effect on the TKE profile. The advanced bubble tracking analysis is completed to post process the simulation results by splitting the bubbles into multiple groups in relation to their distance to the fuel rod surfaces and the distance to the domain inlet. Trends are observed in the population migration, relative velocity and deformability of the bubbles in relation to the complex geometry. Overall, when scaling between high- and low-pressure conditions, the simulations presented in this work display that turbulence closure relations generated for low-pressure flows in subchannel geometries can be applied to both high- and low-pressure conditions. The lowpressure numerical data largely replicated the overall behavior within the high-pressure simulations, specifically when comparing the nearwall behavior of the turbulent in the HP and LP2 cases. However, observing the specifics of bubble dynamics requires more careful practice. The large bubble size at low pressures results in large discrepancies between the population migration and bubble deformation when the geometry of the case is conserved. To ensure the correct bubble behavior is noted when scaling the pressure, scaling the geometry with the size of the bubble is advisable. Specifically, the deformability of the bubbles in the LP2 case replicated the trends of the HP case closely, but the bubbles still had more deformation. The wealth of numerical data that can be collected from these simulations can shed light into the complex two-phase flow behavior within LWR geometries and greatly contribute to the development of closure models with large pressure and temperature application ranges.
References Bolotnov, I.A., 2013. Influence of bubbles on the turbulence anisotropy. J. Fluids Eng. 135 (5), 051301. Fang, J., Bolotnov, I.A., 2017. Bubble tracking analysis of PWR two-phase flow simulations based on the level set method. Nucl. Eng. Design 323, 68–77. Fang, J., Thomas, A.M., Bolotnov, I.A., 2013. In: Transactions of the 2013 ANS Winter Meeting, pp. 1613–1615. Fang, J., Cambareri, J.J., Brown, C.S., Feng, J., Gouws, A., Li, M., Bolotnov, I.A., 2018. Direct numerical simulation of reactor two-phase flows enabled by high-performance computing. Nucl. Eng. Design 330, 409–419. Fang, J., Cambareri, J.J., Rasquin, M., Gouws, A., Balakrishnan, R., Jansen, K.E., Bolotnov, I.A., 2019. Interface tracking investigation of geometric effects on the bubbly flow in PWR subchannels. Nucl. Sci. Eng. 193 (1–2), 46–62. Fletcher, C.D., Schultz, R.R., 1995. RELAP5/MOD3 code manual volume V: User’s Guidelines. Idaho National Engineering Laboratory, Lockheed Idaho Technologies Company, Idaho Falls, Idaho, 83415. Garrett, C., Li, M., Farmer, D., 2000. The connection between bubble size spectra and energy dissipation rates in the upper ocean. J. Phys. Oceanography 30 (9), 2163–2171. Gluck, M., 2008. Validation of the sub-channel code F-COBRA-TF: part I. Recalculation of single-phase and two-phase pressure loss measurements. Nucl. Eng. Des. 238 (9), 2308–2316. Hari, S., Hassan, Y.A., 2002. Improvement of the subcooled boiling model for low-pressure conditions in thermal-hydraulic codes. Nucl. Eng. Design 216 (1–3), 139–152. Ibanez, D.A., Seol, E.S., Smith C.W., Shephard, M.S., 2015. PUMI: parallel unstructured mesh infrastructure. ACM Transactions on Mathematical Software. Jansen, K.E., 1999. A stabilized finite element method for computing turbulence. Comput. Methods Appl. Mech. Eng. 174 (3), 299–317. Jansen, K.E., Johan, Z., Hughes, T.J., 1993. Implementation of a one-equation turbulence model within a stabilized finite element formulation of a symmetric advective-diffusive system. Comput. Methods Appl. Mech. Eng. 105 (3), 405–433. Končar, B., Kljenak, I., Mavko, B., 2004. Modelling of local two-phase flow parameters in upward subcooled flow boiling at low pressure. Int. J. Heat Mass Transfer 47 (6–7), 1499–1513. Kroeger, P.G., Zuber, N., 1968. An analysis of the effects of various parameters on the average void fractions in subcooled boiling. Int. J. Heat Mass Transfer 11 (2), 211–233. Moser, R.D., Kim, J., Mansour, N.N., 1999. Direct numerical simulation of turbulent channel flow up to Re τ= 590. Physics Fluids 11 (4), 943–945. Nagrath, S., Jansen, K.E., Lahey, R.T., 2005. Computation of incompressible bubble dynamics with a stabilized finite element level set method. Comput. Methods Appl. Mech. Eng. 194 (42), 4565–4587. Odar, F., Murray, C., Shumway, R., Bolander, M., Barber, D., Mahaffy, J., 2003. TRAC/ RELAP Advanced Computational Engine (TRACE) V4. 0 User's Manual. US NRC. Rasquin, M., Smith, C., Chitale, K., Seol, E.S., Matthews, B.A., Martin, J.L., Sahni, O., Loy, R.M., Shephard, M.S., Jansen, K.E., 2014. Scalable implicit flow solver for realistic wing simulations with flow control. Comput. Sci. Eng. 16 (6), 13–21. Sussman, M., Fatemi, E., Smereka, P., Osher, S., 1998. An improved level set method for incompressible two-phase flows. Comput. Fluids 27 (5), 663–680. Sussman, M., Almgren, A.S., Bell, J.B., Colella, P., Howell, L.H., Welcome, M.L., 1999. An adaptive level set approach for incompressible two-phase flows. J. Comput. Phys. 148 (1), 81–124. Tu, J.Y., Yeoh, G.H., 2002. On numerical modelling of low-pressure subcooled boiling flows. Int. J. Heat Mass Transfer 45 (6), 1197–1209. Turinsky, P.J., Kothe, D.B., 2016. Modeling and simulation challenges pursued by the Consortium for Advanced Simulation of Light Water Reactors (CASL). J. Comput. Phys. 313, 367–376. Whiting, C.H., Jansen, K.E., 2001. A stabilized finite element method for the incompressible Navier-Stokes equations using a hierarchical basis. Int. J. Numer. Methods Fluids 35, 93–116. Whiting, C.H., Jansen, K.E., 2001. A stabilized finite element method for the incompressible Navier-Stokes equations using a hierarchical basis. Int. J. Numer. Methods Fluids 35 (1), 93–116. Wilcox, D.C., 2002. Turbulence Modeling for CFD, 2nd ed. DCW Industries, La Canada, CA. Woodruff, W.L., Hanan, N.A., Matos, J.E., 1997. A comparison of the RELAP3/MOD3 and PAET/ANL codes with the experimental transient data from the SPERT-IV D-12/25 series. Technical Report of Argonne National Laboratory, Argonne. Zeitoun, O., Shoukri, M., 1997. Axial void fraction profile in low pressure subcooled flow boiling. Int. J. Heat Mass Transfer 40 (4), 869–879.
Declaration of Competing Interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The authors would like to acknowledge the support from the Consortium for Advanced Simulation of Light Water Reactors (CASL) (http://www.casl.gov). CASL is an Energy Innovation Hub for Modeling and Simulation of Nuclear Reactors under U.S. Department of Energy (DOE) [grant number DE-AC05-00OR22725]. The co-author, Dr. Jun Fang, also wants to acknowledge his current employer Argonne National Laboratory for the support related to this work. Argonne is operated by UChicago Argonne, LLC, for the U.S. DOE under contract DE-AC02-06CH11357. This research utilized computational resources provided by Argonne Leadership Computing Facility (ALCF), which is a DOE Office of Science User Facility supported under Contract DE-AC0206CH11357. The solutions presented made use of software components provided by Altair Engineering, Simmetrix (MeshSim and SimModeler),
12
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
Simulation scaling studies of reactor core two-phase flow using direct numerical simulation
T
⁎
Joseph J. Cambareria, , Jun Fangb, Igor A. Bolotnova a b
Department of Nuclear Engineering, North Carolina State University, Raleigh, NC 27695, United States Nuclear Science and Engineering Division, Argonne National Laboratory, Lemont, IL 60439, United States
A R T I C LE I N FO
A B S T R A C T
Keywords: DNS Interface tracking PWR subchannel Scaling studies
Tremendous growth in supercomputing power in recent years has resulted in the emergence of high-resolution flow analysis methods as an advanced research tool to evaluate single and two-phase flow behavior. In particular, unstructured mesh-based methods have been applied to analyze flows in complex reactor core geometries, including those of light water reactors (LWR). The finite-element based code, PHASTA, is utilized to perform large-scale simulations of two-phase bubbly flows in LWR geometries. Given the large computational cost of direct numerical simulation (DNS) coupled with interface tracking methods (ITM), typical domains encompass a portion of a single subchannel. In the presented research, the state-of-the-art analysis of turbulent two-phase flows in complex LWR subchannel geometries are demonstrated at both prototypical reactor parameters as well as scaled low pressure conditions. Three different cases are studied, a high-pressure simulation in prototypical reactor subchannel geometry, a low-pressure case in prototypical geometry and a final low-pressure case in a geometry scaled up to conserve the ratio between the bubble size and the domain pitch. Utilizing advanced statistical processing tools, these simulation conditions are compared to shed light on the relevancy of two-phase flow characteristics given the significant differences between LWR and low-pressure conditions. These findings can lead to the generation of useful guiding principles when researchers need to scale the two-phase flow behavior captured at low pressure and temperature conditions to those at reactor operating conditions.
1. Introduction Due to the high heat transfer coefficients observed within the subcooled boiling regime, the latter is one of the most efficient forms of convective heat transfer. This is in part due to the turbulence intensity generated by the cyclic generation and condensation of vapor bubbles within the flow. However, this regime is limited by the critical heat flux of the system, where past this point highly reduced heat transfer coefficients are observed, resulting in surface burnout or fuel rod damage in nuclear reactor cores. Consideration of the critical heat flux is particularly important in the thermal–hydraulic (TH) and safety analysis of nuclear reactor systems to maintain fuel rod integrity during normal operation as well as accident scenarios. With this concern in mind, the accurate prediction and modeling of turbulent two-phase flows is desired to guide LWR design (Turinsky and Kothe, 2016). This drives the need for expanding the understanding of the physical processes and mechanisms behind turbulent two-phase flows at both high- and lowpressure conditions. For this purpose, multiple TH codes have been developed to support
⁎
the efficient and safe design of nuclear reactors, such as RELAP (Fletcher and Schultz, 1995), TRACE (Odar et al., 2003) and COBRA-TF (Gluck, 2008). While highly cost efficient, the drawback of these TH codes is that they still rely on 1-D or 0-D two-phase closure models. This reliance means that the predictive capabilities of these thermal-hydraulics codes are heavily dependent on the accuracy and applicability of the applied closure relations. Due to the extreme and complex nature of LWR operating parameters, it is highly challenging and expensive to perform high-resolution (interface capturing) experimental measurements at pressures and temperatures observed in nuclear reactors. This results in a lack of quality experimental data for high fidelity model development. For this reason, many subcooled boiling models rely on experimental data from simplified geometries, such as that of Kroeger and Zuber (1968) or the scaling of low-pressure refrigerants to the high pressure and temperature conditions typical of reactor operation. However, this introduces new complexities such as the applicability of scaling relations between low-pressure and high-pressure turbulent bubbly flows and the restrictions of the geometry studied. The breakdown and lack of scalability of these models to low pressure subcooled
Corresponding author. E-mail addresses: [email protected] (J.J. Cambareri), [email protected] (J. Fang), [email protected] (I.A. Bolotnov).
https://doi.org/10.1016/j.nucengdes.2019.110435 Received 30 May 2019; Received in revised form 22 October 2019; Accepted 11 November 2019 0029-5493/ © 2019 Elsevier B.V. All rights reserved.
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
adequately handle turbulence and bubble interactions with and around complex geometries, such as spacer grids or mixing vanes within a LWR subchannel (Fang et al., 2018). Introducing the novel bubble tracking method, which was recently developed, allows for the collection of numerical data on the individual bubble level for the use in advanced statistical post processing (Fang et al., 2013). With the ability to have highly scalable performance on massively parallel systems, up to 786x1024 processors on the IBM Blue Gene/Q Mira system at Argonne National Laboratory, PHASTA is a viable research tool for the study and advanced modeling of turbulent two-phase flows (Rasquin et al., 2014). Mesh generation employs the use of Chef, a derivative of the open source parallel unstructured meshing infrastructure (PUMI) (Ibanez et al., 2015), as well as utilizing tools from Simmetrix, Inc. PHASTA, coupled with Chef meshing tools, is able to work with hexahedral, tetrahedral and mixed finite element meshes along with being able to perform adaptive mesh refinement and coarsening.
boiling conditions for water have also been reported, where unsatisfactory predictions of void fraction distribution are observed. Woodruff et al. (1997) observed this behavior when testing RELAP5 against experimental transient data for the SPERTIV-D-12/25 test. There have also been efforts to improve the modeling of subcooled boiling flows. Hari and Hassan (2002) worked to enhance the predictive capabilities of RELAP5 by implementing the heat partitioning model of Zeitoun and Shoukri (1997) and modifying the interfacial condensation rate. The deficiencies of subcooled boiling models for computational fluids dynamics (CFD) have also been reported by Tu and Yeoh (2002) and Končar et al. (2004), who each proposed multidimensional twofluid subcooled boiling models for the general purpose CFD code, ANSYS-CFX. This indicates a need for further investigation to create guiding principles for researchers when scaling between low-pressure conditions and conditions typical to reactor operation. Generally, better closure relations are required to define the wall heat flux partitioning, the bubble departure diameter and bubble dynamics within the flow. Recently, the expansion of high performance computing (HPC) facilities has led to the viability of utilizing DNS to supplement experimental data when developing accurate two-phase closure models for advanced reactor design (Fang et al., 2018; Fang and Bolotnov, 2017). While DNS is still one of the most computationally expensive methods, the achievable scale and mesh resolution of DNS has observed tremendous growth. Given adequate spatial and temporal resolution, DNS can represent all relevant scales of turbulence without the incorporation of any turbulence models. Coupled with interface tracking methods, DNS can provide high resolution numerical data for turbulent single- and two-phase flows (Moser et al., 1999; Bolotnov, 2013). This ability to collect reliable numerical data has led to DNS becoming an accepted source to supplement experimental data. This work focuses on the study of how the system pressure affects bubble dynamics and interactions with core structures, such as spacer grids and mixing vanes, within LWR geometries. The LWR subchannel geometry is selected to analyze these effects due to the importance of understanding bubble behavior near the mixing vane region and its importance in the calculation of heat transfer mechanisms during normal reactor operation and accident scenarios. Three different simulations are considered. The first recreates the conditions that are expected for typical reactor operation within an LWR subchannel geometry and the analysis of this case has been previously presented in literature (Fang et al., 2018). The second and third simulations scale the case parameters to those expected at atmospheric conditions, with one simulation conserving the geometry and the other scaling the geometry to match the new bubble diameter to hydraulic diameter ratio. In both simulations, the bubble deformability, through the Eötvös number (Eo), are preserved by increasing the bubble diameter with the differing fluid properties. For each simulation, a previously developed bubble tracking method was applied to track the statistics of individual bubbles within the flow (Fang et al., 2013). These new simulations and state-of-the-art statistical analysis of turbulent two-phase flows at both low-pressure conditions and prototypical reactor operating conditions will help to shed light onto relevant bubbly flow behavior in LWR geometries and the impact of spacer grid and mixing vanes on two-phase models and closure laws.
2.1. Flow solver The conditions considered within this research describe an isothermal and incompressible system. PHASTA solves the weak form of the Incompressible Navier-Stokes (INS) equations in three dimensions through a stabilized FEM. The spatial and temporal discretization for the INS equations within PHASTA have been previously discussed at length by Whiting and Jansen (2001) as well as Nagrath et al. (2005). The governing equations for continuity and momentum in the strong form of the INS equations are given by:
∂ui =0 ∂x i
ρ
∂τij ∂p ∂ui ∂u + fi + + ρuj i = − ∂x j ∂x i ∂t ∂x j
(1)
(2)
here, ui denotes the velocity for the i-th dimension with i = 1, 2 and 3, ρ the fluid density, p the static pressure, τij as the viscous stress tensor and finally fi represents the body force in the i-th coordinate, including the gravitational and surface tension forces. This final term includes both the gravitational and surface tension forces. The viscous stress tensor for an incompressible flow of a Newtonian fluid is connected to the dynamic viscosity, μ , of the fluid and the strain rate tensor, Sij .
∂uj ⎞ ∂u τij = 2μSij = μ ⎜⎛ i + ⎟ ∂x i ⎠ ∂ x j ⎝
(3)
Expanding from single- to two-phase capabilities, PHASTA employs the level-set method to distinguish between the liquid and gas phase (Sussman et al., 1998, 1999). The level-set method creates a precise, continuous, signed distance function, φ , which represents the bubble interface at the zero of the level-set distance field. Thus, φ = 0 represents the interfacial boundary between the phases. By convention, a positive value of the level-set distance field indicates the liquid phase, while the gas phase is denoted by a negative value. The scalar, φ , is advected with the fluid according to the advection equation.
∂φ ∂φ =0 + ui ∂x i ∂t
2. Numerical methods
(4)
To account for the discontinuity created by the fluid properties at the interface, a smoothed Heaviside kernel function, H I µ̂ , is introduced (Sussman et al., 1999). This function assumes a finite interfacial thickness, and the properties over the interface, the density and viscosity of each phase, are smoothed to remove the jump discontinuity, as presented below.
The finite element method (FEM) research code, PHASTA (Jansen, 1999), was used to perform the direct numerical simulations. PHASTA is a parallel, hierarchic, higher-order accurate, adaptive, stabilized transient analysis flow solver for the solution to both incompressible and compressible flows. The flow solver has been shown to be a highly effective research tool for a broad range of turbulent flow conditions and scales, having been applied to RANS, LES, Detached Eddy Simulation (DES) and DNS cases (Jansen et al., 1993; Whiting and Jansen, 2001). With its ability to process unstructured grids, PHASTA can 2
ρ (φ) = ρ1 Hε (φ) + ρ2 (1 − Hε (φ))
(5)
μ (φ) = μ1 Hε (φ) + μ 2 (1 − Hε (φ))
(6)
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
2.2. Bubble tracking method Throughout the history of DNS, there have been multiple techniques for numerical data collection and analysis. Single- and two-phase quantities had been collected through the implementation of static virtual probes arrayed in 2-D planes within the computational domain. The virtual probes planes are placed normal to the streamwise direction, with probe locations set at specific distances from the wall. For each distance to the wall value there is a specific number of virtual probes where similar statistical behavior is expected to occur. Over the course of the simulation, these probes collect instantaneous numerical data, such as the phase velocities or void fraction. Statistical averaging in time and space is then used to compute two-phase turbulent flow parameters, such as the mean flow velocity, void fraction distributions, Reynolds stresses, turbulent kinetic energy or the turbulence dissipation rate of the liquid phase. While the virtual probe statistical analysis approach has made major contributions to the development and study of turbulent closure models (Bolotnov, 2013), especially for singlephase turbulent flows, the efficient and useful numerical data collection for two-phase turbulent flows presents more of a challenge. The virtual probe arrays are in static planar locations, removing their ability to adequately capture the spatial movement patterns and behavior of bubbles in two-phase flows, especially in complex geometries such as those with spacer grids and mixing vanes. For this reason, the development and implementation of a bubble tracking method was sought after. Developed by Fang et al. (2013), a scalar marker field is initialized along with the level set distance field. This scalar field individually tracks and collects the instantaneous numerical data of each separate bubble within the computational domain. Each bubble is assigned a unique ID, which is advected with the level set to dynamically track the bubbles simulated. Additionally, numerical data from the surrounding liquid shell is collected to assess the local flow conditions, such as the liquid shear or relative velocity, which the bubble is encountering. The bubble information that can be extracted from this method is highly valuable for the generation and improvement of models and closure laws, such as the drag or lift coefficient, for use in lower resolution computational fluid dynamics (CFD) to better predict turbulent bubble flows.
a)
b) Fig. 1. a) Computational domain setup for the three cases (front walls hidden to display inner structures) and b) location of the measurement planes within the domain.
spacer grid and mixing vanes, as displayed in Fig. 1 below. Here, the front walls of the domain are hidden to display the inner structures of the spacer grid and mixing vanes. The length of the subchannel domain is 4.43 rod diameters in length (40.5 mm). This domain length was selected to balance the quality of numerical statistics and the computational cost of the simulations. Two measurement planes of static virtual probes are located within the computational domain, set at locations of 0.86 and 2.28 rod diameters from the end of the mixing vanes (2.74 and 4.16 rod diameters from the domain origin), as displayed by the red dashed lines in Fig. 1. Three different simulation setups are considered in this investigation: a 1% void fraction (262 bubbles) highpressure case previously presented in literature (Fang et al., 2018, 2019), a simulation that conserves the prototypical subchannel geometry while scaling to low-pressure conditions and a final simulation which scales the geometry size to conserve the bubble diameter to pitch ratio for low-pressure conditions. As a note, for the following comparison the simulations will be referred to as HP for the high-pressure simulation, LP1 for the case conserving the initial geometry and scaling to low-pressure conditions and LP2 for the simulation which scales both the pressure and geometric conditions. Periodic boundary conditions are applied at the inlet and outlet of the domain as well as the transverse faces to allow for the turbulence and bubble interactions with adjacent subchannels without the need for meshing those regions. Along the surface of the fuel rods, spacer grid and mixing vanes the no-slip wall boundary condition is applied.
3. Simulation setup When considering the thermal–hydraulic analysis of an LWR, the subchannel geometry plays a critical role in the overall reactor core design. As coolant flows through the subchannel, it will carry heat generated from the reactor fuel rods to the steam generators. During prototypical PWR operating conditions, subcooled boiling can occur in the upper region of the reactor core, resulting in the production of steam bubbles within the subchannel coolant. The inclusion of vapor bubbles within the flow typically results in an increase of heat transfer efficiency because of the bubble induced turbulence and the subcooled boiling phenomenon. On top of this, the inclusion of complex geometries within the subchannel, such as spacer grids and mixing vanes, can greatly affect the solution of the coolant flow field, enhancing the turbulent intensity of the coolant and increasing the heat transfer and heat removal from the fuel rod surfaces. While both vapor production and complex geometries further complicate the heat transfer mechanisms at play, the full impact these effects on two-phase flows within LWR subchannels are not yet fully understood due to the difficulty of obtaining high-resolution experimental measurements at the extreme temperatures and pressures observed in the reactor core environment.
3.2. Parameter scaling to low pressure simulations Detailed results of the 1% void fraction case at conditions approaching prototypical reactor operating conditions have been presented in literature (Fang et al., 2018, 2019). For this simulation, the fluid properties were set by using the saturated properties of water and vapor at 573 K. The saturation temperature selected for this work is associate with a corresponding system pressure of 8.58 MPa, displayed in Table 1, which is equivalent to conditions observed within an LWR. However, in the push to create guiding principles for researchers to follow when scaling between pressure conditions, the following research focuses on the comparison between the numerical data gathered from those investigations to DNS of turbulent bubble flows completed at low-pressure conditions. For the low-pressure simulations, the system pressure was selected to be 0.1 MPa and the fluid properties are based
3.1. Computational domain design Given the available computational resources, and knowing the importance of the subchannel geometry, the computational domain generated for this work includes a single subchannel with a reduced size 3
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
Table 1 Property comparison at prototypical reactor conditions and low-pressure conditions. Parameter/Case
High-pressure
Low-pressure
Unit
Pressure Saturation Temperature Fluid Density (liquid; vapor) Fluid Dynamic Viscosity (liquid; vapor) Surface Tension Coefficient Gravitational Acceleration Bubble Diameter Eötvös Number
8.58 573 712; 46 8.59E−05; 1.97E−05 0.014 9.81 0.65 0.25
0.10 373 959; 0.59 2.82E−04; 1.23E−05 0.059 9.81 1.25 0.25
MPa K kg/m3 Pa s N/m m/s2 mm –
Table 2 Overview of significant geometric and computational properties. Parameter/Case
HP
LPC1
LPC2
Bubble diameter (mm) Hinze Scale (mm) Rod diameter (mm) Rod pitch (mm) Domain height (mm) Mesh cell size in bulk region (μm) Mesh cell size of the first boundary layer (μm) Total number of mesh cells Number of bubbles resolved Mean velocity (m/s) Flow Reynolds number
0.65 0.82 9.14 12.60 40.50 29.30 8.14
1.25 1.57 9.14 12.60 40.50 57.21 7.15
1.25 2.04 17.66 24.33 78.22 61.68 7.71
1.10 billion 262 0.75 8.1E04
0.29 billion 37 0.75 3.3E04
1.13 billion 262 0.75 6.4E04
Fig. 3. Scaled virtual probe locations for the LP2 simulation.
of the bubbles had to be approximately doubled to ensure the same value of the Eötvös number. By guaranteeing the conservation of this dimensionless characteristic, the behavior of the bubbles in the prototypical case will be most readily compared to those in the low-pressure simulations.
on the corresponding saturated properties of water and vapor. Of main interest in this study was the comparison of bubble behavior and dynamics around complex LWR geometries (e.g. spacer grids and mixing vanes) when scaling to low pressure conditions. For this reason, the bubble deformability, through the Eötvös number, presented in (7) below, was conserved. Here, the Eötvös number is calculated using the equivalent diameter, db , of a perfectly spherical bubble. From the previous studies of turbulent bubbly flows in LWR subchannels, the resulting average Eötvös number of the resolved bubbles was 0.25. Because of the large change in simulation properties resulting from the change in pressure, mainly due to the influence of the larger surface tension coefficient at atmospheric pressures, the diameter
a)
Eo =
Δρgdb2 σ
(7)
It should be noted that the Reynolds number of the original case was not conserved throughout each simulation. At higher Reynolds numbers the bubbles are more likely to breakup as the influence of the external forces upon the bubble overcome the surface tension forces. The impact of this can be generally classified with the calculation of the Hinze scale, given in (8) below (Garrett et al., 2000). Here, c is a model constant with a value of 0.363 and ε is the dissipation rate. If the bubble radius is greater than the Hinze scale it is more likely to breakup under the flow conditions. For the high- and low-pressure simulations the value of the Hinze scale are listed in Table 2 below. In each of these
b)
Fig. 2. Initial bubble size and geometry comparison of the high- and low-pressure cases for a) the conserved geometry and b) the scaled geometry. 4
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b) Fig. 4. Snapshots for the LP1, HP and LP2 simulations (left to right) at a) the two-phase initialization (0 BDT) and b) 30 BDT with the bubbles delineated by their IDs and the slices displaying the instantaneous velocity.
the turbulence of the flow. Since the both the high- and low-pressure bubble radius falls below the Hinze scale, the lack of Reynolds number scaling is not anticipated to have a large impact on the bubble
simulations the bubble size is much larger than the Kolmogorov length scale (on the order of 10 µm for these simulations) and smaller than the Hinze scale, so minimal breakup of the bubbles is to be expected due to 5
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
Fig. 5. The instantaneous liquid velocity field at the first measurement plane with the location of the outer tip of the mixing vanes marked (dashed black line).
a)
b)
Fig. 6. The liquid a) mean streamwise velocity and b) velocity tangential to the fuel rod surfaces at the first measurement plane location.
a)
b)
Fig. 7. The turbulent kinetic energy profile of the liquid phase versus the mixing vane tip location (black dashed line) for the a) first measurement plane and b) second measurement plane.
bubbles to only 37. To create a more direct comparison of the highpressure case, another low-pressure simulation was initialized, but the geometric properties of the high-pressure case were scaled to conserve the bubble diameter to hydraulic diameter ratio of the original geometry. To conserve this ratio, the new domain was scaled linearly by a factor of 1.93 times the size of the original domain. Additionally, to replicate the original high-pressure simulation conditions as close as possible, the initial bubble locations of the 1% void fraction case were scaled with the same geometric scaling factor. A comparison of
dynamics.
aH = c (σ / ρ)3/5ε −2/5
(8)
With the large change in the bubble diameter resulting from the conservation of the Eötvös number, the ratio of the bubble diameter to the geometric conditions of the original computational domain were greatly affected. To maintain the 1% void fraction initialization that was applied to the high-pressure simulation, the number of bubbles in the atmospheric pressure case had to be greatly reduced from 262 6
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b)
c) Fig. 8. Anisotropy distribution comparison for the a) HP and LP1 case, b) HP and LP2 case and c) the difference between the data sets where a one-to-one fit (dashed black line) is presented with lines for ± 20% (dotted black lines).
Now, the smaller turbulent eddies can be developed on the finer mesh for a more time efficient process. The refinement process can continue until the required turbulent scales are achieved. For each simulation, the desired turbulent profile and flow rate is accomplished by adjusting the flow driving pressure gradient in the streamwise direction until it balances the wall drag on the flow to achieve the desired flow rate. The inclusion of the spacer grid and mixing vanes creates a highly efficient process for the development of single-phase turbulence. The mean flow velocity is set at 0.75 m/s for each void fraction case, which corresponds to a hydraulic Reynolds number of 80,737 for the HP case, 33,097 for LP1 and 63,812 for LP2. The gravitational force is applied in the direction opposite to the flow to create the conditions of vertical upward channel flow. Once the single-phase turbulence reaches a statistically steady state, the two-phase conditions are initialized. Fig. 4 exhibits the initial bubble locations, colored by their ID numbers, for the HP simulation and the low-pressure cases, as well as displaying the initial instantons velocity field. Fig. 4 also displays a further two-phase time step for the simulations at approximately 30 bubble diameter time (BDT) based on the HP case. The BDT is a nondimensional unit that allows for a more direct comparison between the cases and is defined as the average time interval that a bubble required to move one diameter in distance. For the presented analysis the BDT of the LP1 and LP2 cases was converted to use the bubble size of the HP case. To ensure the quality of the statistics collected, the two-phase numerical data was collected once a quasi-steady state was reached for each simulation. High quality statistics are a necessity for numerical simulations and one of the easiest methods to accomplish this is to
important geometric and computational properties for the three cases are displayed in Table 2 below. To better visualize the change in the geometric and simulation parameters, Fig. 2 shows the high and low pressure bubble sizes with both the conserved and scaled geometries.
4. Results and discussion All of the cases were run on the IBM Blue Gene/Q Mira supercomputer at Argonne National Laboratory (ANL). The new LP1 and LP2 cases were carried out on a significantly larger computational allocation than previously awarded, so the mesh quality, based on the number of elements across a bubble diameter, was increased to better resolve the interfacial behavior. From the previous 1% void fraction case, the number of elements across a bubble diameter was increased from 21 to approximately 24. Both the conventional static array of probes and the novel bubble tracking algorithm data collection methods were applied for numerical data collection in the simulation scaling studies. Two measurement planes of static virtual probes are placed within the computational domain, set at locations of 0.86 and 2.28 rod diameters from the end of the mixing vanes for each case. For the LP2 case, the locations of the probes were scaled with the same geometric ratio as before. The resulting probe plane is presented in Fig. 3. The numerical data collected by the probe arrays is supplemented with that from the bubble tracking algorithm. For an efficient large-scale PHASTA simulation workflow, the large turbulent eddies are initially developed on a coarse mesh before the mesh is uniformly refined and the solution migrated onto a finer mesh. 7
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b)
c) Fig. 9. Anisotropy distribution of the off-diagonal components for the a) HP and LP1 case, b) HP and LP2 case and c) the difference between the data sets where a one-to-one fit (dashed black line) is presented with lines for ± 20% (dotted black lines). Table 3 Comparison of group bounds for the high- and low-pressure simulations. HP
LP1
LP2
Group Number
Spatial range (by wall distance, unit: mm)
Group Number
Spatial range (by wall distance, unit: mm)
Group Number
Spatial range (by wall distance, unit: mm)
F1 F2 F3 F4 F5 F6
(0.0, 0.51) (0.51, 0.98) (0.98, 1.41) (1.41, 1.82) (1.82, 2.47) (2.47, 4.34)
FR1
(0.0, 0.98)
FR2
(0.98, 1.82)
FR3
(1.82, 4.34)
F1 F2 F3 F4 F5 F6
(0.0, 0.99) (0.99, 1.89) (1.89, 2.71) (2.71, 3.51) (3.51, 4.77) (4.77, 8.69)
increase the sample size of the numerical data set. Due to the large computational cost and the transient nature of turbulent DNS simulations, this is not always an option though. To balance the quality of statistics gathered and the computational cost, special attention is given to the convergent behavior of the time averaged statistics. Three metrics were used to assess this condition; the convergence of the mean streamwise velocity, the balance of the bubble population in relation to the distance to the wall and the development of the bubble relative velocity. The numerical data collected from the virtual probe planes is averaged across multiple time windows, and the convergence of the streamwise velocity profile is evaluated. The second two methods employ the bubble tracking numerical data. The bubble population distribution and relative velocity can be dynamically tracked over the course of the entire simulation. The solution is assumed to reach a quasi-steady state condition when the relative vapor/liquid velocity has stabilized in the streamwise direction and the distribution of bubbles in
Fig. 10. The bubble tracking groups bounds for HP and LP2 (solid lines) and LP1 (dotted bands) in a quarter subchannel cross-section.
relation to the distance to the wall has stabilized. For this work, it was identified that approximately 40 BDT was required to reach this quasisteady state, at which point the streamwise velocity profile had also converged. Due to the fuel rod geometry, the local wall normal coordinates at the virtual probe locations change in reference to the global coordinate system. This presents a unique challenge in processing the numerical data, where the wall distance at two opposite fuel rod surfaces has a 8
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b)
Fig. 11. Evolution of the bubble population for the a) near wall and b) channel center distance to the fuel rod surface groups as a function of the distance from the domain inlet for the HP (solid lines), LP1 (dotted lines) and LP2 (dashed lines) cases.
from subchannels without the spacer grid and mixing vane do not show this trend, where a flatter profile is seen outside of the boundary layer region (Fang and Bolotnov, 2017). By the new coordinate system convention, the tangential direction is positive in a counterclockwise direction along the surface of the fuel rods. Due to the dominant clockwise rotation of the coolant, the tangential velocity shows a large negative velocity between the locations of 0.4 and 0.8 r/R. Within this region the mixing vanes have the largest influence upon the fluid flow. On top of the mean velocity profile, the static virtual probe planes also can process the velocity fluctuations in the domain to calculate the turbulent kinetic energy (TKE) profile for the liquid phase. The impact of the internal structures is also prominent on the turbulent kinetic energy profiles. Compared to a simple subchannel without an internal structures (Fang and Bolotnov, 2017), the presence of the spacer grid and mixing vanes significantly increases the magnitude of the TKE profile in the coolant flow, as seen in Fig. 7. Within the center of the subchannel, as well as in the near wall region for LP1, the low-pressure cases observe higher values for the TKE. However, for both of the low-
zero value. To account for this aspect of the subchannel geometry, the Cartesian coordinates within each quarter of the subchannel are adjusted into a streamwise, tangential and wall normal component system. To produce the time averaged statistics, the instantaneous numerical data from the virtual probes located at each distance to the wall value are averaged over a specific time window. For this study, the numerical data collected by the virtual probes within the rod gaps was excluded to reduce the influence of the cross flow on the statistics. With the inclusion of internal structures in the subchannel, especially resulting from the mixing vane geometry, the instantaneous liquid velocity field reveals a centrifugal pattern as observed in Fig. 5. The mixing vanes produce the rotating flow in the core of the subchannel, forcefully pushing the flow in the direction of the outer regions. This has the additional effect of generating a lateral flow at the surface of the fuel rods, enhancing cross flow and mixing between subchannels. The mean streamwise velocity, depicted in Fig. 6, also displays this effect, where a reduction in the liquid velocity in the center of the subchannel is observed. Previous mean velocity profiles 9
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b)
Fig. 12. The bubble relative velocity for the a) near wall and b) channel center distance to the fuel rod surface groups as a function of the distance from the domain inlet for the HP (solid lines), LP1 (dotted lines) and LP2 (dashed lines) cases.
prototypical reactor and low-pressure conditions using state-of-the-art statistical analysis based on the dynamic bubble tracking method. Through the implementation of the bubble tracking algorithm, detailed bubble information is collected over the entire lifetime of a bubble during the two-phase flow simulations. The collected bubble tracking data is divided into several bubble groups based on their distance to the fuel rod surface and the distance to the domain inlet for a more detailed study of the bubble dynamics involved. For this study of the HP case, the bubble numerical data was split into six groups based on their distance to the wall and fourteen groups based on the distance to the domain inlet. The group outer bounds were set such that each group has an equivalent cross sectional area, and the group bounds are presented in Table 3. Adjustments were also made for the group bounds in the LP1 case due to the size increase of the bubbles within the conserved subchannel geometry. Specifically, the number of distance to the fuel rod groups was reduced to three due to the size of the bubbles in relation to the geometry. A visual representation of the group bounds within a quarter of the subchannel cross section is displayed in Fig. 10 below. Within this figure, the solid lines with the dotted bands represent where the group bounds for the HP and LP2 cases coincide with that of the LP1 case. Since the initial bubble size is larger than that of the viscous sublayer for all cases, limited numerical bubble data can be gathered for the near wall relations. However, important bubble behavior and dynamics can still be evaluated within the turbulent log-law region of the subchannel geometry. To track the specific bubbles that occupy each group, the centroid of the bubble is recorded, and its location is compared in relation to the group bounds. While this is a fairly simple approach, future advancements of this technique can account for when
pressure cases the influence of the complex geometry fades away faster than that of the high-pressure case. Future studies will investigate the impact of bubble induced turbulence to isolate the dominate effect changing the TKE profile. With the application of DNS in this work, each term within the Reynolds stress tensor can be computed for the liquid phase. Specifically, the directional components of the Reynolds stress tensor −
are calculated for the streamwise (u'u' ), normal to the fuel rod surface − ' '
−
(v v ) and tangential to the fuel rod surface (w 'w ' ) components. The twophase flow anisotropy distribution can be computed by normalizing these values with the turbulent kinetic energy. For a flat plat boundary layer the ratio of the normal stresses is typically found to be −
−
−
u'u': v 'v ': w 'w ' = 8: 4: 6, so the subchannel anisotropic values are scaled by a factor of 9 for an easier comparison (Wilcox, 2002), however with the inclusion of the complex geometry in the subchannel domain this value is expected to differ. Fig. 8 presents the comparison of the highand low-pressure data as well as the difference between the values. The comparison in Fig. 8 presents all of the directional components of the Reynolds stress tensor in relation to a one-to-one and a difference of ± 20%. Overall, the LP2 case appears to more closely match the turbulent anisotropy distribution of the HP case. The off-diagonal components of the Reynolds stress tensor also play an important role in the overall characterization of the subchannel turbulence. Fig. 9 presents the comparison of the off-diagonal Reynolds stress tensor terms for the HP, LP1 and LP2 cases. While a wealth of information can be calculated from the turbulent statistics of DNS, this work mainly focuses on the comparison of 10
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
a)
b)
Fig. 13. The bubble deformability for the a) near wall and b) channel center distance to the fuel rod surface groups as a function of the distance from the domain inlet for the HP (solid lines), LP1 (dotted lines) and LP2 (dashed lines) cases.
Further analysis that can be completed with the bubble tracking method involve the tracking of novel statistics such as the bubble deformability factor. The bubble deformability factor, γd , is defined as the ratio of the minimum value of the level set field found within a bubble to the radius of a sphere with the equivalent volume as the considered bubble. In short, this term quantifies how deformed a bubble is, with an ideal sphere having a value of γd = 1.0 . Fig. 13 displays the effect that the complex geometry have upon the shape of bubbles within the flow, where bubbles entering and leaving the spacer grid and mixing vane region, especially in the core of the flow and close to the wall, are notably more deformed than those outside of the complex geometry region of influence. Within the spacer grid region this deformation is likely caused by flow confinement, and after the mixing vanes the bubbles tend to deform under the influence of the turbulent flow. Both of these effects decay as the bubbles move downstream from the region of influence of the complex geometry. For both low-pressure cases, the bubbles are significantly more deformed than that of the high-pressure case, with LP1 observing the highest amounts of deformation. This deformability is most likely a result of the relatively large bubble size moving though the complex geometry coupled with the increased weight of the low-pressure fluid. Additionally, for the central bubble groups, the bubbles tend to return to a more spherical shape after leaving the area of influence generated by the spacer grid and mixing vanes within the channel.
a bubble is moving across a group bound and indicate the bubbles that are entering and leaving a group. Separating the bubbles into groups based on the distance to the fuel rod allows for the evaluation of the population distribution and bubble relative velocity throughout the subchannel over the course of the bubbly flow simulations, as seen in Fig. 11 and Fig. 12. For this work the relative velocity convention is taken as the bubble velocity minus the liquid velocity. Within these plots the solid vertical black lines indicate the beginning and the end of the spacer grid region, while the dashed vertical black line specifies the end of the mixing vanes. This method allows the observation of this migration trend throughout the length of the subchannel and the evolution of bubble dynamics as they are affected by the spacer grid and mixing vanes. Generally, upward two-phase flows typically have a void fraction profile that peaks closer to the wall, and previous subchannel simulations excluding the spacer grid and mixing vanes demonstrate this migration trend of bubbles towards the wall (Fang and Bolotnov, 2017). However, the centrifugal effects of the mixing vanes alter the dynamics of the flow, forcing bubbles to migrate into the center of the subchannel after they encounter the mixing vanes. Notably, after the mixing vane region, the most central groups (F6 and FR3) see large increase in bubble population due to the turbulent mixing effects of the mixing vanes. This increase in population is accompanied by a sharp decrease in the relative velocity of the bubbles found within those groups, although this decrease is observed for every distance to the fuel rod surface group. Both the LP1 and LP2 cases observe a larger migration of bubbles into the core of the flow, however a smaller relative velocity drop when compared to the HP case.
5. Conclusion In the presented research, turbulent bubbly flows within LWR 11
Nuclear Engineering and Design 358 (2020) 110435
J.J. Cambareri, et al.
SCOREC at RPI (Chef) and Kitware (ParaView).
subchannels are investigated at conditions prototypical to reactor operation and those found in low-pressure scenarios. Using high-fidelity, state-of-the-art interface tracking simulations, all relevant scales of liquid turbulence are resolved, and the interface is captured with the level-set method. Three different cases are considered to carry out the scaling study of the effects of system pressure upon bubble dynamics within the flow. Using static virtual planes, the impact of the complex geometry can be observed on both the mean flow velocities and the TKE profile. Further investigations into the effect of bubble induced turbulence and the change in liquid properties are planned to examine its effect on the TKE profile. The advanced bubble tracking analysis is completed to post process the simulation results by splitting the bubbles into multiple groups in relation to their distance to the fuel rod surfaces and the distance to the domain inlet. Trends are observed in the population migration, relative velocity and deformability of the bubbles in relation to the complex geometry. Overall, when scaling between high- and low-pressure conditions, the simulations presented in this work display that turbulence closure relations generated for low-pressure flows in subchannel geometries can be applied to both high- and low-pressure conditions. The lowpressure numerical data largely replicated the overall behavior within the high-pressure simulations, specifically when comparing the nearwall behavior of the turbulent in the HP and LP2 cases. However, observing the specifics of bubble dynamics requires more careful practice. The large bubble size at low pressures results in large discrepancies between the population migration and bubble deformation when the geometry of the case is conserved. To ensure the correct bubble behavior is noted when scaling the pressure, scaling the geometry with the size of the bubble is advisable. Specifically, the deformability of the bubbles in the LP2 case replicated the trends of the HP case closely, but the bubbles still had more deformation. The wealth of numerical data that can be collected from these simulations can shed light into the complex two-phase flow behavior within LWR geometries and greatly contribute to the development of closure models with large pressure and temperature application ranges.
References Bolotnov, I.A., 2013. Influence of bubbles on the turbulence anisotropy. J. Fluids Eng. 135 (5), 051301. Fang, J., Bolotnov, I.A., 2017. Bubble tracking analysis of PWR two-phase flow simulations based on the level set method. Nucl. Eng. Design 323, 68–77. Fang, J., Thomas, A.M., Bolotnov, I.A., 2013. In: Transactions of the 2013 ANS Winter Meeting, pp. 1613–1615. Fang, J., Cambareri, J.J., Brown, C.S., Feng, J., Gouws, A., Li, M., Bolotnov, I.A., 2018. Direct numerical simulation of reactor two-phase flows enabled by high-performance computing. Nucl. Eng. Design 330, 409–419. Fang, J., Cambareri, J.J., Rasquin, M., Gouws, A., Balakrishnan, R., Jansen, K.E., Bolotnov, I.A., 2019. Interface tracking investigation of geometric effects on the bubbly flow in PWR subchannels. Nucl. Sci. Eng. 193 (1–2), 46–62. Fletcher, C.D., Schultz, R.R., 1995. RELAP5/MOD3 code manual volume V: User’s Guidelines. Idaho National Engineering Laboratory, Lockheed Idaho Technologies Company, Idaho Falls, Idaho, 83415. Garrett, C., Li, M., Farmer, D., 2000. The connection between bubble size spectra and energy dissipation rates in the upper ocean. J. Phys. Oceanography 30 (9), 2163–2171. Gluck, M., 2008. Validation of the sub-channel code F-COBRA-TF: part I. Recalculation of single-phase and two-phase pressure loss measurements. Nucl. Eng. Des. 238 (9), 2308–2316. Hari, S., Hassan, Y.A., 2002. Improvement of the subcooled boiling model for low-pressure conditions in thermal-hydraulic codes. Nucl. Eng. Design 216 (1–3), 139–152. Ibanez, D.A., Seol, E.S., Smith C.W., Shephard, M.S., 2015. PUMI: parallel unstructured mesh infrastructure. ACM Transactions on Mathematical Software. Jansen, K.E., 1999. A stabilized finite element method for computing turbulence. Comput. Methods Appl. Mech. Eng. 174 (3), 299–317. Jansen, K.E., Johan, Z., Hughes, T.J., 1993. Implementation of a one-equation turbulence model within a stabilized finite element formulation of a symmetric advective-diffusive system. Comput. Methods Appl. Mech. Eng. 105 (3), 405–433. Končar, B., Kljenak, I., Mavko, B., 2004. Modelling of local two-phase flow parameters in upward subcooled flow boiling at low pressure. Int. J. Heat Mass Transfer 47 (6–7), 1499–1513. Kroeger, P.G., Zuber, N., 1968. An analysis of the effects of various parameters on the average void fractions in subcooled boiling. Int. J. Heat Mass Transfer 11 (2), 211–233. Moser, R.D., Kim, J., Mansour, N.N., 1999. Direct numerical simulation of turbulent channel flow up to Re τ= 590. Physics Fluids 11 (4), 943–945. Nagrath, S., Jansen, K.E., Lahey, R.T., 2005. Computation of incompressible bubble dynamics with a stabilized finite element level set method. Comput. Methods Appl. Mech. Eng. 194 (42), 4565–4587. Odar, F., Murray, C., Shumway, R., Bolander, M., Barber, D., Mahaffy, J., 2003. TRAC/ RELAP Advanced Computational Engine (TRACE) V4. 0 User's Manual. US NRC. Rasquin, M., Smith, C., Chitale, K., Seol, E.S., Matthews, B.A., Martin, J.L., Sahni, O., Loy, R.M., Shephard, M.S., Jansen, K.E., 2014. Scalable implicit flow solver for realistic wing simulations with flow control. Comput. Sci. Eng. 16 (6), 13–21. Sussman, M., Fatemi, E., Smereka, P., Osher, S., 1998. An improved level set method for incompressible two-phase flows. Comput. Fluids 27 (5), 663–680. Sussman, M., Almgren, A.S., Bell, J.B., Colella, P., Howell, L.H., Welcome, M.L., 1999. An adaptive level set approach for incompressible two-phase flows. J. Comput. Phys. 148 (1), 81–124. Tu, J.Y., Yeoh, G.H., 2002. On numerical modelling of low-pressure subcooled boiling flows. Int. J. Heat Mass Transfer 45 (6), 1197–1209. Turinsky, P.J., Kothe, D.B., 2016. Modeling and simulation challenges pursued by the Consortium for Advanced Simulation of Light Water Reactors (CASL). J. Comput. Phys. 313, 367–376. Whiting, C.H., Jansen, K.E., 2001. A stabilized finite element method for the incompressible Navier-Stokes equations using a hierarchical basis. Int. J. Numer. Methods Fluids 35, 93–116. Whiting, C.H., Jansen, K.E., 2001. A stabilized finite element method for the incompressible Navier-Stokes equations using a hierarchical basis. Int. J. Numer. Methods Fluids 35 (1), 93–116. Wilcox, D.C., 2002. Turbulence Modeling for CFD, 2nd ed. DCW Industries, La Canada, CA. Woodruff, W.L., Hanan, N.A., Matos, J.E., 1997. A comparison of the RELAP3/MOD3 and PAET/ANL codes with the experimental transient data from the SPERT-IV D-12/25 series. Technical Report of Argonne National Laboratory, Argonne. Zeitoun, O., Shoukri, M., 1997. Axial void fraction profile in low pressure subcooled flow boiling. Int. J. Heat Mass Transfer 40 (4), 869–879.
Declaration of Competing Interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The authors would like to acknowledge the support from the Consortium for Advanced Simulation of Light Water Reactors (CASL) (http://www.casl.gov). CASL is an Energy Innovation Hub for Modeling and Simulation of Nuclear Reactors under U.S. Department of Energy (DOE) [grant number DE-AC05-00OR22725]. The co-author, Dr. Jun Fang, also wants to acknowledge his current employer Argonne National Laboratory for the support related to this work. Argonne is operated by UChicago Argonne, LLC, for the U.S. DOE under contract DE-AC02-06CH11357. This research utilized computational resources provided by Argonne Leadership Computing Facility (ALCF), which is a DOE Office of Science User Facility supported under Contract DE-AC0206CH11357. The solutions presented made use of software components provided by Altair Engineering, Simmetrix (MeshSim and SimModeler),
12