Issue 
EPJ Nuclear Sci. Technol.
Volume 6, 2020
Euratom Research and Training in 2019: the Awards collection



Article Number  6  
Number of page(s)  7  
DOI  https://doi.org/10.1051/epjn/2019057  
Published online  07 February 2020 
https://doi.org/10.1051/epjn/2019057
Regular Article
Hydraulic and statistical study of metastable phenomena in PWR rod bundles
CEA Saclay, DEN/DANS/DM2S/STMF/LMSF, 91191 GifsurYvette, France
^{*} email: florian.muller99@gmail.com
Received:
30
October
2019
Accepted:
15
November
2019
Published online: 7 February 2020
The analysis of fuel rod bundle flows constitutes a key element of PressurizedWater Reactors (PWR) safety studies. The present work aims at improving our understanding of nefarious reorganisation phenomena observed by numerous studies in the flow largescale structures. 3D simulations allowed identifying two distinct reorganisations consisting in a sign change for either a transverse velocity in rodtorod gaps or for a subchannel vortex. A Taylor “frozen turbulence” hypothesis was adopted to model the evolution of largescale 3D structures as transported2D. A statistical method was applied to the 2D field to determine its thermodynamically stable states through an optimization problem. Similarities were obtained between the PWR coherent structures and the stable states in a simplified 2D geometry. Further, 2D simulations allowed identifying two possible flow bifurcations, each related to one of the reorganisations observed in 3D simulations, laying the foundations for a physical explanation of this phenomenon.
© F. Muller, published by EDP Sciences, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Insufficient flow thermal mixing in the rod bundles within a PressurizedWater Reactor (PWR) can lead to a boiling crisis, which is nefarious for the reactor operational safety. Mixing grids are typically used to enhance the thermal mixing inside fuel arrays, mostly through the intensification of the secondary flows. These secondary flows have a tendency to organize into largescale structures in the plane normal to the tube axes. Numerous experimental or numerical studies have shown the existence of reorganisation phenomena in the transverse flow largescale structures (see the review in appendix A from Kang and Hassan [1]). In particular, the AGATE experimental results [2] featured a global 90° rotation of the crossflow pattern between the near and the far wake of a mixing grid. This reorganisation is shown in Figure 1: a sketch of a 5 × 5 rod bundle fitted with a mixing grid, as well as colour plots of the transverse flow downstream a mixing grid are shown. The crossflow is aligned with a 45° angle in the first one, but has rotated by 90° in the second one.
Such reorganisations are very significant for reactor safety studies: due to the points with zero crossflow they induce, they lead to drops in the thermal mixing as demonstrated by Shen et al. [3], which can pose a serious risk to the PWR reactor operational safety. This work aimed at improving our understanding of these phenomena, both for enhancing their characterization and for identifying their origin, with the longterm goal of developing smallscales models suited for this type of flow.
Little concrete information can be found in the literature on the reorganisation phenomenon. This is due to the lack of highquality experiment where the phenomenon typically occurs, and to the fact that, among the variety of ComputationalFluid Dynamics (CFD) numerical simulations performed for rod bundle flows, few were conducted for the entire rod bundle axial span and with highfidelity turbulence models. Attempting a new method, we focused on a physical approach to the reorganisation phenomenon by proposing an original method for the study of the rod bundle flow. A similarity was noticed between PWR rod bundle flow reorganisations and some phenomena typically experienced by quasi2D geophysical flows, such as the Jupiter Red Spot or the Gulf Stream [4]. Indeed, the latter sometimes display important changes in their structure, leading to oscillations between very distinct solutions. These phenomena can be interpreted as phase transitions between different equilibrium states which become metastable.
In order to study the reorganisation phenomena from the perspective of this similarity with 2D flows, the following steps were taken. 3D simulations were first performed in order to decompose the largescale reorganisations into local inversions, and to justify a Taylor “frozen turbulence” hypothesis, as described in Section 2. Section 3 details the 2D statistical method that was applied in simplified geometries with obstacles based on this hypothesis. The stable states obtained through this method are then used in 2D free decay simulations in Section 4, highlighting similarities between their phase transitions and the 3D reorganisation phenomenon. Section 5 provides a conclusion on the physicality of the reorganisation phenomenon.
Fig. 1 Sketch of a rod bundle and colour plots of the transverse velocity field downstream the mixing grid. 
2 3D CFD simulations
We performed 3D CFD simulations in rod bundles using the TrioCFD [5] code along with a WALE subgrid scale turbulence model, in a tetrahedral mesh of 34.4 million cells and using in parallel 1000 CPU cores. Details on the numerical performance of this code were discussed by Angeli et al. [6]. The flow was characterized by a Reynolds number Re = 96 000. Crossflow reorganisations were observed in the simulation results, as depicted in Figure 2: the averaged transverse flow at 30 mm downstream the mixing grid in the first plot displays a 45° orientation, while that at 150 mm in the second plot is aligned with the 135° diagonal. A significant impact of the numerical schemes used was also noted, which is a topic still under investigation by the community. Most importantly, we could distinguish two particular types of local crossflow reorganisations.
The crossflow can undergo an inversion of the transverse velocity between two rods, or a sign change for a vortex in the centre of a subchannel (the vacant space between four rods). Different combinations of “velocity inversions” around a subchannel (e.g. one to four inversions) can lead to a variety of pattern changes, one of which is a global 90° rotation of the crossflow pattern. This type of inversion was already noted in previous work by Shen et al. [3].
Besides, we also observed in the 3D simulation results multiple subchannels featuring a crossflow pattern change without any gap flow inversion. Instead, the circulation inside these subchannels seemed to vary greatly. This increase consisted either in a sign change for the subchannel vortex if one was present, or the apparition of a vortex in a pure shear flow. One of the possible largescale pattern changes obtained was the 90° rotation of the subchannel crossflow pattern.
Remarkably, the same final patterns can be obtained after each of these two reorganisations through two different flow evolutions, providing a more nuanced view of the possible reorganisations in rod bundle crossflows.
In order to tackle the reorganisation phenomenon, it was decided to adopt an approach based on a statistical fluid mechanics point of view and relying on a Taylor “frozen turbulence” hypothesis. This hypothesis is commonly used in meteorological flows (see Higgins et al. [7]) where local turbulent eddies are advected by a mean wind and are considered “frozen” at a fixed position. Under this hypothesis, one can consider that a rod bundle crossflow behaves as a quasi2D flow transported by the axial component of the velocity field. This hypothesis allows decomposing the 3D velocity field u(x, y, z, t) into:
Here is a uniform axial component, and u ′ (x, y, z, t) denotes the 3D turbulent fluctuations. The (x, y) components form a purely advected transverse 2D part that is assumed to abide by the 2D Navier–Stokes equations. In order to verify this hypothesis, the flow must display very small variations of the axial component compared to that of the transverse flow field , and a high scale difference between the axial and transverse components.
These criteria were mostly fulfilled outside of boundary layers and past the near wake of the mixing grid, notably with a scale ratio between the axial and transverse components of up to 30. Although deviations from the hypothesis are observed and should be investigated further, such a decomposition allowed us to apply statistical methods inspired by geophysical flow studies to the transverse 2D part.
Fig. 2 Averaged transverse velocity field from a LES simulation of the flow in a 5 × 5 rods fuel assembly, at 30 mm (left) and 150 mm (right) downstream a mixing grid. 
3 2D statistical approach
Past studies on 2D turbulent flows through the lens of statistical fluid mechanics have tackled a broad set of physical fields, from thin oceanic layers to the Red Spot in the Jovian atmosphere and experimental soap films. A crucial part of the study of such 2D flows is the prediction of the steady stable structures emerging from given initial characteristics and depending on the domain geometry. In the framework of PWR rod bundle flows, this prediction consists in the determination of stable patterns for the transverse flow, following “initial conditions” set by its passage through the mixing grid.
Instead of an exhaustive application of a dynamical stability criterion to the infinite set of steady solutions to the 2D Euler equations, an equivalence can be used between such a dynamical stability criterion and a constrained optimization problem, usually on a measure of the entropy with conservation of a varying set of Euler invariants, following the work from Miller et al. [8,9]. This allows for the use of several tools coming from mathematical optimization theory, and to directly calculate the expected steady state in a given 2D physical configuration.
Multiple theories have been devised to determine these stable states, leading to various constrained optimization problems. We have focused on the MinimumEnstrophyPrinciple (MEP) due to its relative simplicity and the linear equations it leads to. This principle was proposed from phenomenological considerations (see Bretherton and Haidvogel [10]), and it states that the stable flow regimes minimize the enstrophy of the system under constant energy and circulation. The enstrophy is defined as , with ω the flow vorticity field and Ω the physical domain.
A resolution method for this problem was devised and proposed by Naso et al. [11], applicable to simplyconnected domain geometries. We adapted this method to geometrical domains with internal obstacles in order to allow the consideration of typical rod bundle crosssection geometries. Each obstacle adds the conservation of the flow circulation around itself and thus another degree of freedom to the problem.
The minimization of the enstrophy under multiple constraints leads to a variational problem solved in our method through a projection on the Laplace eigenbasis. Scanning the parameter space and interpolating for constant energy, circulation and circulation around the obstacles allows determining the stable solutions for any combination of the integral constraints. Details on the derivation of the solutions and on the method followed to explore the solution ensemble are available in Muller et al. [12].
Our method has been validated first through numerical convergence studies and recovery of literature results [13], before being applied in the case of a square domain with two diagonally opposed obstacles. This simple geometry was designed as a basic model for a PWR rod bundle crosssection. In order to reduce the number of degrees of freedom in the problem, we assumed equal circulations around the two obstacles, justified by the typical crossflow patterns observed downstream mixing grids in PWR rod bundle flows and shown in Figure 1.
This approach resulted in the stability map shown in Figure 3. Therein, Λ^{2} = Γ^{2}/(2E) with E the total flow energy combines the energy and circulation into one parameter, while Γ_{1} is the circulation around the obstacles in this geometry. The thick black line indicates the minimalenstrophy state, solution of our multiplyconstrained optimization problem. Under the MEP, this solution is a stable flow for given values of the parameters(Γ, E, Γ_{1}). The light and dark grey areas, respectively, indicate areas of lowerenstrophy for 1 and 2vortices solutions, but their enstrophy is systematically higher than that of the solutions on the black line. In the graph, β is related to the energy E and a _{1} to the boundary condition on the obstacles. Interestingly, a single solution is stable for, while two branches are stable for. Further, the single branch features a central vortex aligned with the diagonal without obstacles and one vortex around each obstacle, while the upper branch for displays a single vortex around both obstacles on the opposite diagonal.
Therefore, two stable regimes differing by a 90° angle were obtained, which recalls the pattern reorganisations observed in PWR rod bundle crossflows. Although in the 2D statistical approach, the parameters (Γ, E, Γ_{1}) are fixed due to their conservation by the Euler equations, they can be dissipated by viscosity in physical systems, leading to changes in the minimalenstrophy solution and thus stable states. Phase transitions can thus occur between these states, which were investigated further in our study.
Fig. 3 Stability map obtained in a square geometry with two obstacles, following the application of the MinimumEnstrophy Principle. Depending on the value of the parameter Λ^{2}, one or two flow regimes can be stable in the domain. 
4 2D freedecay CFD simulations
To further validate the methodology and the computational code for the statistical approach and check if a link could be established between the statistical approach and the CFD calculations, the case of the free relaxation of turbulent 2D flows has been studied afterwards. The objective was to compare the flow reached in the longrun of a 2D CFD simulation and the stable state predicted by the statistical approach, from initial conditions either random or designed to investigate particular flow regimes.
The simulations were set up to feature as little diffusion as possible in order to conserve energy relatively well when compared to the spontaneous enstrophy dissipation. Similarly to the 3D simulations, the 2D simulations conducted in this study were based on the code TrioCFD [5]. They did not use any turbulence model. In order to allow the capture of most flow scales, moderate Reynolds numbers were chosen (Re = 3000). When choosing the numerical schemes, we opted as much as possible for the less dissipative ones: only slightly upstream convection schemes, a thirdorder RungeKutta time advancement method and a resolution of the pressure equation through a Cholesky method [6]. Regarding the boundary conditions, we used “stressfree” boundary conditions, which in practice amounted to a nullity condition on the normal velocity component at the boundary. We found this method to be the closest one to the physical framework of the statistical approach, where the existence of a nonzero circulation Γ_{1} around the internal obstacles requires in turn the tangential velocity at the boundary of the obstacles to be nonzero.
We first validated the process by recovering final states in square and rectangular domain geometries in accordance with statistical mechanics literature results, and by capturing plausible turbulent cascades for the kinetic energy. The decay process obtained in a square geometry is shown in Figure 4. In this geometry, the MinimumEnstrophy Principle approach (see Chavanis et al. [13]) leads to a single central vortex as the only stable solution. Having initiated our simulation with an array of TaylorGreen vortices, it underwent a phase of strong mixing, before converging towards the expected singlevortex stable state.
In the case of the square domain with two obstacles, the simulations were tailored to have initial conditions spanning the (Λ^{2}, Γ_{1}/Γ) parameter space. This was performed through the initialization of the simulations with carefully designed analytical functions. Let (O, e_{r}, e_{θ})define a polar coordinate system with O the square center as origin point, e_{r} a unitary vector along the position vector r and e_{θ} a unitary vector orthogonal to r. The analytical functions for the simulation initial conditions were built through the addition of the following terms:

a rigid body rotation term u = Ωre_{θ} (with Ω an arbitrary angular velocity and r the distance to the origin O) leading to a background vortex in the entire square domain;

a single vortex around each obstacle, built for each one from the combination of a 2 × 2 array of TaylorGreen vortices with damping functions depending on the distance from the obstacle centre.
The rotation sign of the obstacle vortices were set opposite to that of the central vortex, and the relative intensity of the two counterrotating phenomena was varied between the different tests in order to adjust the initial value of Γ _{1}/Γ. The simulations were performed using two mesh refinements, both with 21 000 or 237 000 triangular cells. These limited resolutions were imposed by the limited resources available due to the large number of cases tested and the long simulated time required in each case. The starting points were all placed in the righthand side of the stability map shown in Figure 3, i.e. in the zone where the system should only present one stable solution. We thus hoped to observe the flow select either of the two branches when following the free decay process and having a decreasing Λ^{2}.
The energy, total circulation and circulation around the obstacles were calculated in order to observe the evolution of the freelydecaying 2D simulations in the (Λ^{2}, Γ_{1}/Γ) stability map. The resulting evolutions are superimposed on the stability map in Figure 5. Note that the flow remained quite symmetrical during the simulations, to the extent that the circulations around the two obstacles remained about equal. This allowed a comparison with the similar case of the square with two obstacles of equal circulations in Section 3.
Among the 30 simulations tested from various initial points, two types of bifurcations were observed, respectively displayed as continuous and dashed lines in Figure 5.

Either the central vortex of negative vorticity was conserved while the vortices around the rods switched their rotation direction, implying a better conservation of Γ than Γ_{1}. This bifurcation is displayed in the velocity vector fields in Figure 6, and is annotated from “1a” to “1d” in Figure 5. Such a bifurcation implies a sign change for the velocity between obstacles and walls, recalling the “velocity inversion” phenomenon observed in our simulations in Section 2 or from the work of Shen et al. [3].

Or the central vortex switched its rotation direction, “forced” by the vortices around the obstacles which in this case are conserved. This bifurcation is displayed in the velocity vector fields in Figure 7 (with a nearzero vorticity for the initial central vortex), and is annotated from “2a” to “2d” in Figure 5. This bifurcation is comparable to the sign change for a subchannel vortex observed in our simulation results for rod bundle flows.
To assess the “pertinence” of the adopted approach based on 2D equations, different approaches were considered to establish links between flows obtained from 2D simulations and the crossflows observed in planes orthogonal to the main axial direction in 3D simulation results.
First, initial and intermediary states from 2D freedecay simulations were inserted as a crossflow into a 3D flow with a large advection velocity in order to observe the decay process in a quasi2D setup. Similarities were observed between the transported2D and 2D evolutions; in particular a comparable final state was observed in some cases. However, this approach encountered significant challenges such as 3D diffusion, computational domain design and boundary layer effects. Second, a 2D simulation with an initial condition as the crossflow taken directly after a 3 × 3 rod bundle mixing grid in a 3D flow has been realized. Parallels could be drawn between its resulting decay process and that observed in the 3D crossflow, but the 2D simulation was observed to reach in the long run a symmetrical state definitely not achievable by the 3D flow due to the shortness of the rod bundle axial span.
The conclusion from these attempts is that the results obtained in 2D are difficult to directly transpose in the full 3D flow. Among others, boundary layer effects and departures from the Taylor hypothesis require further investigation before clear quantitative parallels can be drawn between the 3D simulation and experimental results and the 2D theoretical and numerical results.
Fig. 4 Colour plot of the vorticity field along a free decay 2D simulation in a square, from an array of TaylorGreen vortices to a minimalenstrophy single vortex. 
Fig. 5 Decay processes (red and blue plots) in a square geometry with two obstacles, within the related stability map (Fig. 3). Two possible bifurcation paths were identified. 
Fig. 6 Decay process in a square geometry with two obstacles and with a type 1 inversion; each vortex around an obstacle switches its sign. The left and right plots correspond, respectively, to points “1b” and “1d” in Figure 5. 
Fig. 7 Decay process in a square geometry with two obstacles and with a type 2 inversion; the central vortex switches its sign. The left and right plots correspond respectively to points “2b” and “2c” in Figure 5. 
5 Conclusion
The various new learnings from this work can be synthesized as the following: the global reorganisation phenomenon observed in PWR rod bundle experiments and numerical simulations can be decomposed at the local level into sign changes for either the rod gap velocity or the subchannel vortex. Decomposing the flow into axial and transverse parts based on a Taylor “frozen turbulence” hypothesis allows studying the transverse flow using tools from 2D statistical fluid mechanics. A resolution method was designed to allow for the identification of stable states in domains with internal obstacles based on the minimization of the system enstrophy. In a square with two obstacles modelling a PWR rod bundle crosssection, two particular flow regimes differing by a 90° rotation were observed, comparably to the PWR crossflow patterns. Freelydecaying 2D simulations in this geometry lead to the identification of two bifurcation mechanisms, each related to an inversion phenomenon in PWR rod bundle flows. This comparison paves the way for a better understanding of reorganisations in nuclear cores coolant flows.
From an industrial point of view, this work should encourage members of the nuclear research community to intensify the investigation of PWR rod bundle crossflow reorganisation phenomena. The phenomena can be observed in most rod bundle experimental results; the disturbances in thermal mixing they entail can raise significant safety issues. Further experimental data should therefore be acquired and published in the open literature. The numerical simulations performed in this work in 2D and 3D should be pursued and refined, in particular in order to improve the links between the 2D and 3D approaches through the inclusion of boundary layer effects and deviation from the quasi2D framework.
As mentioned in Section 4, results from the 2D analysis are challenging to directly transpose into a 3D framework. However, the idea is that a mixed approach should be considered. In the long term, a hybrid 2D–3D model could provide firstestimate results for 3D flows in a fraction of the computational costs, by combining:

a 3D simulation of the flow in the near wake of a mixing grid, where the Taylor “frozen turbulence” hypothesis is invalid;

a 2D simulation in the rest of the rod bundle using as initial condition the outlet crossflow of the 3D simulation.
Numerous roadblocks still need to be lifted before such a model can be applied into industrial problems, but it could provide interesting perspectives in the study of PWR rod bundle flows.
References
 S.K. Kang, Y.A. Hassan, Computational fluid dynamics (CFD) round robin benchmark for a pressurized water reactor (PWR) rod bundle, Nucl. Eng. Des. 301, 204 (2016) [CrossRef] [Google Scholar]
 U. Bieder, F. Falk, G. Fauchet, LES analysis of the flow in a simplified PWR assembly with mixing grid, Prog. Nucl. Energy 75, 15 (2014) [CrossRef] [Google Scholar]
 Y. Fen Shen, Z. Dong Cao, Q. Gang Lu, An investigation of crossflow mixing effect caused by grid spacer with mixing blades in a rod bundle, Nucl. Eng. Des. 125, 111 (1991) [CrossRef] [Google Scholar]
 M.J. Schmeits, H.A. Dijkstra, Bimodal behavior of the Kuroshio and the Gulf Stream, J. Phys. Oceanogr. 31, 3435 (2001) [CrossRef] [Google Scholar]
 TrioCFD: http://wwwtriou.cea.fr/ [Google Scholar]
 P.E. Angeli, M.A. Puscas, G. Fauchet, A. Cartalade, FVCA8 Benchmark for the Stokes and NavierStokes Equations with the TrioCFD Code—Benchmark Session, in Finite Volumes for Complex Applications VIII − Methods and Theoretical Aspects , edited by C. Cancès, P. Omnes (Springer International Publishing, Cham, 2017), p. 181 [CrossRef] [Google Scholar]
 C.W. Higgins, M. Froidevaux, V. Simeonov, N. Vercauteren, C. Barry, M.B. Parlange, The effect of scale on the applicability of Taylor's frozen turbulence hypothesis in the atmospheric boundary layer, Boundary Layer Meteorol. 143 , 379 (2012) [CrossRef] [Google Scholar]
 J. Miller, Statistical mechanics of Euler equations in two dimensions. Phys. Rev. Lett. 65, 2137 (1990) [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 R. Robert, J. Sommeria, Statistical equilibrium states for twodimensional flows, J. Fluid Mech. 229, 291 (1991) [CrossRef] [MathSciNet] [Google Scholar]
 F.P. Bretherton, D.B. Haidvogel, Twodimensional turbulence above topography, J. Fluid Mech. 78 , 129 (1976) [CrossRef] [Google Scholar]
 A. Naso, P.H. Chavanis, B. Dubrulle, Statistical mechanics of twodimensional Euler flows and minimum enstrophy states, Eur. Phys. J. B 77, 187 (2010) [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 F. Muller, A. Burbeau, B.J. Gréa, P. Sagaut, Minimum enstrophy principle for twodimensional inviscid flows around obstacles, Phys. Rev. E 99 , 023105 (2019) [CrossRef] [Google Scholar]
 P.H. Chavanis, J. Sommeria, Classification of selforganized vortices in twodimensional turbulence: the case of a bounded domain, J. Fluid Mech. 314, 267 (1996) [NASA ADS] [CrossRef] [Google Scholar]
Cite this article as: Florian Muller, Hydraulic and statistical study of metastable phenomena in PWR rod bundles, EPJ Nuclear Sci. Technol. 6, 6 (2020)
All Figures
Fig. 1 Sketch of a rod bundle and colour plots of the transverse velocity field downstream the mixing grid. 

In the text 
Fig. 2 Averaged transverse velocity field from a LES simulation of the flow in a 5 × 5 rods fuel assembly, at 30 mm (left) and 150 mm (right) downstream a mixing grid. 

In the text 
Fig. 3 Stability map obtained in a square geometry with two obstacles, following the application of the MinimumEnstrophy Principle. Depending on the value of the parameter Λ^{2}, one or two flow regimes can be stable in the domain. 

In the text 
Fig. 4 Colour plot of the vorticity field along a free decay 2D simulation in a square, from an array of TaylorGreen vortices to a minimalenstrophy single vortex. 

In the text 
Fig. 5 Decay processes (red and blue plots) in a square geometry with two obstacles, within the related stability map (Fig. 3). Two possible bifurcation paths were identified. 

In the text 
Fig. 6 Decay process in a square geometry with two obstacles and with a type 1 inversion; each vortex around an obstacle switches its sign. The left and right plots correspond, respectively, to points “1b” and “1d” in Figure 5. 

In the text 
Fig. 7 Decay process in a square geometry with two obstacles and with a type 2 inversion; the central vortex switches its sign. The left and right plots correspond respectively to points “2b” and “2c” in Figure 5. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.