Nội dung

EPJ Nuclear Sci. Technol. 6, 6 (2020) © F. Muller, published by EDP Sciences, 2020 https://doi.org/10.1051/epjn/2019057 Nuclear Sciences & Technologies Available online at: https://www.epj-n.org REGULAR ARTICLE Hydraulic and statistical study of metastable phenomena in PWR rod bundles Florian Muller* CEA Saclay, DEN/DANS/DM2S/STMF/LMSF, 91191 Gif-sur-Yvette, France Received: 30 October 2019 / Accepted: 15 November 2019 Abstract. The analysis of fuel rod bundle flows constitutes a key element of Pressurized-Water Reactors (PWR) safety studies. The present work aims at improving our understanding of nefarious reorganisation phenomena observed by numerous studies in the flow large-scale structures. 3D simulations allowed identifying two distinct reorganisations consisting in a sign change for either a transverse velocity in rod-to-rod gaps or for a subchannel vortex. A Taylor “frozen turbulence” hypothesis was adopted to model the evolution of large-scale 3D structures as transported-2D. 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. 1 Introduction Insufficient flow thermal mixing in the rod bundles within a Pressurized-Water 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 large-scale 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 large-scale 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 cross-flow 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 cross-flow 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 cross-flow 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, * e-mail: florian.muller99@gmail.com both for enhancing their characterization and for identifying their origin, with the long-term goal of developing small-scales 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 high-quality experiment where the phenomenon typically occurs, and to the fact that, among the variety of Computational-Fluid Dynamics (CFD) numerical simulations performed for rod bundle flows, few were conducted for the entire rod bundle axial span and with high-fidelity 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 quasi-2D 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 large-scale 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 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. 2 F. Muller: EPJ Nuclear Sci. Technol. 6, 6 (2020) 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. 2 3D CFD simulations We performed 3D CFD simulations in rod bundles using the TrioCFD [5] code along with a WALE sub-grid scale Fig. 1. Sketch of a rod bundle and colour plots of the transverse velocity field downstream the mixing grid. 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. Cross-flow 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 cross-flow reorganisations. The cross-flow 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 cross-flow 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 sub-channels featuring a cross-flow pattern change without any gap flow inversion. Instead, the circulation inside these sub-channels 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 large-scale pattern changes obtained was the 90° rotation of the subchannel cross-flow 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 cross-flows. 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. F. Muller: EPJ Nuclear Sci. Technol. 6, 6 (2020) 3 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 cross-flow behaves as a quasi-2D flow transported by the axial component of the velocity field. This hypothesis allows decomposing the 3D velocity field u(x, y, z, t) into: uðx; y; z; tÞ ¼ uz z þ ux;y ðx; y; uz tÞ þ u0 ðx; y; z; tÞ: Here uz is a uniform axial component, and u0 (x, y, z, t) denotes the 3D turbulent fluctuations. The (x, y) components form a purely advected transverse 2D part ux;y ðx; y; uz tÞ 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 uz compared to that of the transverse flow field ux;y , 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. 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 Minimum-EnstrophyPrinciple (MEP) due to its relative simplicity and the linear equations it leads to. This principle was proposed Fig. 3. Stability map obtained in a square geometry with two obstacles, following the application of the Minimum-Enstrophy Principle. Depending on the value of the parameter L2, one or two flow regimes can be stable in the domain. 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 G2 ¼ ∫v2 dr, with v the flow vorticity field and V the V 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 cross-section 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 cross-flow 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, L2 = G2/(2E) with E the total flow energy combines the energy and circulation into one 4 F. Muller: EPJ Nuclear Sci. Technol. 6, 6 (2020) Fig. 4. Colour plot of the vorticity field along a free decay 2D simulation in a square, from an array of Taylor-Green vortices to a minimal-enstrophy single vortex. parameter, while G1 is the circulation around the obstacles in this geometry. The thick black line indicates the minimal-enstrophy state, solution of our multiplyconstrained optimization problem. Under the MEP, this solution is a stable flow for given values of the parameters (G, E, G1). The light and dark grey areas, respectively, indicate areas of lower-enstrophy for 1- and 2-vortices solutions, but their enstrophy is systematically higher than that of the solutions on the black line. In the graph, b is related to the energy E and a1 to the boundary condition on the obstacles. Interestingly, a single solution is stable forL2 > L2crit , while two branches are stable forL2 < L2crit . 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 L2 < L2crit 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 cross-flows. Although in the 2D statistical approach, the parameters (G, E, G1) are fixed due to their conservation by the Euler equations, they can be dissipated by viscosity in physical systems, leading to changes in the minimal-enstrophy solution and thus stable states. Phase transitions can thus occur between these states, which were investigated further in our study. 4 2D free-decay 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 long-run 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 third-order Runge-Kutta time advancement method and a resolution of the pressure equation through a Cholesky method [6]. Regarding the boundary conditions, we used “stress-free” 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 non-zero circulation G1 around the internal obstacles requires in turn the tangential velocity at the boundary of the obstacles to be non-zero. 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 single-vortex stable state. In the case of the square domain with two obstacles, the simulations were tailored to have initial conditions spanning the (L2, G1/G) parameter space. This was performed through the initialization of the simulations with carefully designed analytical functions. Let (O, er, eu) define a polar coordinate system with O the square center as origin point, er a unitary vector along the position vector r and eu 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 = Vreu (with V 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 Taylor-Green vortices with damping functions depending on the distance from the obstacle centre. F. Muller: EPJ Nuclear Sci. Technol. 6, 6 (2020) The rotation sign of the obstacle vortices were set opposite to that of the central vortex, and the relative intensity of the two counter-rotating phenomena was varied between the different tests in order to adjust the initial value of G 1/G. 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 right-hand 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. 5 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 L2. The energy, total circulation and circulation around the obstacles were calculated in order to observe the evolution of the freely-decaying 2D simulations in the (L2, G1/G) 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 G than G1. 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 near-zero vorticity for the initial central vortex), and is annotated 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. 6 F. Muller: EPJ Nuclear Sci. Technol. 6, 6 (2020) 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. 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 cross-flows 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 cross-flow into a 3D flow with a large advection velocity in order to observe the decay process in a quasi-2D setup. Similarities were observed between the transported-2D 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 cross-flow 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 cross-flow, 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. 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 cross-section, two particular flow regimes differing by a 90° rotation were observed, comparably to the PWR cross-flow patterns. Freely-decaying 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 cross-flow 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 F. Muller: EPJ Nuclear Sci. Technol. 6, 6 (2020) 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 quasi-2D 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 first-estimate 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 cross-flow 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 1. 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) 2. 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) 7 3. 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) 4. M.J. Schmeits, H.A. Dijkstra, Bimodal behavior of the Kuroshio and the Gulf Stream, J. Phys. Oceanogr. 31, 3435 (2001) 5. TrioCFD: http://www-trio-u.cea.fr/ 6. P.-E. Angeli, M.-A. Puscas, G. Fauchet, A. Cartalade, FVCA8 Benchmark for the Stokes and Navier-Stokes 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 7. 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) 8. J. Miller, Statistical mechanics of Euler equations in two dimensions. Phys. Rev. Lett. 65, 2137 (1990) 9. R. Robert, J. Sommeria, Statistical equilibrium states for two-dimensional flows, J. Fluid Mech. 229, 291 (1991) 10. F.P. Bretherton, D.B. Haidvogel, Two-dimensional turbulence above topography, J. Fluid Mech. 78, 129 (1976) 11. A. Naso, P.H. Chavanis, B. Dubrulle, Statistical mechanics of two-dimensional Euler flows and minimum enstrophy states, Eur. Phys. J. B 77, 187 (2010) 12. F. Muller, A. Burbeau, B.-J. Gréa, P. Sagaut, Minimum enstrophy principle for two-dimensional inviscid flows around obstacles, Phys. Rev. E 99, 023105 (2019) 13. P.H. Chavanis, J. Sommeria, Classification of self-organized vortices in two-dimensional turbulence: the case of a bounded domain, J. Fluid Mech. 314, 267 (1996) Cite this article as: Florian Muller, Hydraulic and statistical study of metastable phenomena in PWR rod bundles, EPJ Nuclear Sci. Technol. 6, 6 (2020)