A coarse-mesh methodology for modelling of single-phase thermal-hydraulics of ESFR innovative assembly design

Nuclear Engineering and Design 355 (2019) 110291

Contents lists available at ScienceDirect

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

A coarse-mesh methodology for modelling of single-phase thermalhydraulics of ESFR innovative assembly design

T

Stefan Radmana, , Carlo Fiorinaa, Konstantin Mikityukb, Andreas Pautza,b ⁎

a b

École Polytechnique Fédérale de Lausanne, Laboratory for Reactor Physics and Systems Behaviour, 1015 Lausanne, Switzerland Paul Scherrer Institut, Nuclear Energy and Safety Department, 5232 Villigen, Switzerland

ARTICLE INFO

ABSTRACT

Keywords: Sodium Fast Reactor SFR Coarse-mesh Porous medium Inter-assembly gap Assembly wrapper Wrapper windows OpenFOAM GeN-Foam

The GeN-Foam multi-physics code enables coupled thermal-hydraulic, neutronic, and thermal-mechanical analysis of fast reactor cores with 3-D unstructured deformable meshes. Advances within the Horizon2020 ESFRSMART project, as well as trends towards higher-fidelity core-wise simulations motivated further developments of the code thermal-hydraulics, based on a coarse-mesh porous medium approach. These developments enable the single-phase characterization of the thermal-hydraulic impact of the inter-assembly gap and assembly wrapper windows in steady state and transient scenarios. The inter-assembly gap constitutes a non-negligible heat sink in accident scenarios in most SFR designs. Conversely, wrapper windows are an ESFR innovative feature meant to alleviate the adverse effects of sodium boiling. However, these will alter the single-phase thermal-hydraulics as well. This work presents development and verification efforts aimed to enable the coarsemesh thermal-hydraulic analysis of these features. An analysis of the impact of these features is performed on a representative bundle of 7 ESFR assemblies.

1. Introduction

developed by the Laboratory for Reactor Physics and Systems Behaviour (LRS) at the EPFL, Switzerland, and the Paul Scherrer Institut (PSI), Switzerland. The code was developed for the 3-D multi-physics analysis of thermal-hydraulics, neutronics and thermal-mechanics at a core level with general unstructured deformable meshes. In order to achieve the required level of modelling flexibility without incurring in excessive computational burdens, the solver thermal-hydraulics is built upon a porous medium based coarse-mesh approach. The ultimate goal is to enable the analysis of reactor cores at an assembly-wise level of detail, yet without a high-fidelity simulation of the intra-assembly phenomenology. It is thus possible to perform calculations on mediumend workstations without the need for supercomputing clusters. This approach allows for the investigation of 3-D effects while limiting the required computational power compared to standard CFD approaches. Phenomena that can hardly be captured by 1-D codes can be thus investigated, such as thermal stratification in pool-type reactor plena (Wahnon et al., 2018), cross-flow in the inter-assembly gap and between assemblies, and natural circulation in loss-of-flow scenarios. A porous medium approach for nuclear applications is not, however, a novelty in its own regard, as testified by various codes such as the legacy BACCHUS-3D/TP (Bottoni et al., 1987) or recent codes such as TWOPORFLO (Chavez et al., 2018), and CUPID (Jeong et al., 2010). What does represent an element of novelty is the implementation of this

Sodium-cooled Fast Reactors (SFRs) are among the nearest term deployable fast reactor technologies explored in the frame of the Generation-IV International Forum (GIF) (OECD-NEA, 2014). This is due to significant operating experience acquired internationally over the past 50 years, although interest in SFR technology varied significantly over said time scale. On a European level, the 90s saw a slowdown in SFR research and investment (Tenchine, 2010), culminating in the shutdown of the Superphenix reactor (CEA, 1978), and the termination of the European Fast Reactor (EFR) project (Lefvre et al., 1996). Eventually, the interest in designing and constructing a novel SFR in Europe was reaffirmed by collaborative projects. These efforts, stemming from the original EFR project, evolved into the Collaborative Project on European Sodium Fast Reactor (CP ESFR) (Fiorini and Vasile, 2011) and the currently ongoing Horizon 2020 project European Sodium Fast Reactor Safety Measures Assessment and Research Tools (ESFR-SMART) (Mikityuk et al., 2017). It is in this context that the present work is set. The analysis of fast reactor systems along with the general trend towards higher fidelity simulations calls for adequate multi-physics approaches. An answer to such need is represented by the OpenFOAM® (OpenCFD Ltd., 2019; Weller and Tabor, 1998; Jasak, 2009) based GeN-Foam multi-physics code (Fiorina et al., 2015), jointly ⁎

Corresponding author. E-mail address: [email protected] (S. Radman).

https://doi.org/10.1016/j.nucengdes.2019.110291 Received 25 March 2019; Received in revised form 29 May 2019; Accepted 11 August 2019 Available online 19 September 2019 0029-5493/ © 2019 Elsevier B.V. All rights reserved.

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

methodology in a modern programming paradigm-based library such as OpenFOAM, which significantly simplifies multi-physics coupling, and its use over general unstructured 3-D meshes. The various physics of the code and their coupling underwent considerable verification efforts against established codes such as TRACE, PARCS, Serpent, MCNP (Fiorina et al., 2015; Fiorina et al., 2016; Fiorina et al., 2017) and analytical solutions (de Oliveira and Mikityuk, 2018) for a variety of cases. Regular sanity checks (Radman et al., 2018) of the code multiphysics capabilities are performed to aid the development process. Recent developments in SFR core design, such as the trend towards a low void reactivity feedback (Bortot et al., 2015), result in a tighter balance between different reactivity feedbacks during accident scenarios. This strengthens the need for higher fidelity simulation tools to further assess the impact of particular SFR design features and their behaviour during transients. One of the innovative features proposed for SFR cores consists in assembly wrapper windows (Mikityuk et al., 2011). These are intended to limit the negative consequences that arise from sodium boiling. With respect to thermal reactors, boiling in fast reactors can lead to a positive reactivity insertion (Sun et al., 2011; Tommasi et al., 2010). Sodium boiling presents an additional challenge as it can lead to a significant flow reduction and exacerbate a ULOF scenario (Boure et al., 1973; Kottowski and Savatteri, 1984). The flow excursion is due to the considerable increase of the friction factor in two-phase flow conditions. This in turn facilitates the downward axial propagation of the sodium void, potentially reaching core regions characterized by a local high positive void reactivity coefficient. A boiling stabilization technique can thus consist in promoting sodium vapour presence only to the upper core region, which can be achieved e.g. with the introduction of openings in the assembly wrappers, which will be referred to as assembly windows. Inter-assembly gap flow is a standard feature in SFRs because of the use of a wrapper tube for every assembly. It can constitute a non-negligible heat sink during partial assembly blockage scenarios, given the excellent thermal properties of liquid sodium (> 60 W/m/K of thermal conductivity) and the favourable heat transfer properties of the assembly steel wrapper. While the role of the inter-assembly gap used to be neglected or strongly approximated in the past, the need for a more detailed representation of the inter-assembly flow has recently been recognized in some investigations (Yue et al., 2018; Droin et al., 2017; Hu and Yu, 2016; Radman et al., 2017) and is being enabled by further developments in 2-D and 3-D modelling capabilities of nuclear codes. The thermal-hydraulic solver of GeN-Foam was expanded to enable the prediction of assembly-gap heat transfer phenomena across the assembly wrapper as well as localized pressure drops at the assembly windows within the context of a coarse-mesh methodology. The thermal-hydraulic impact of said features was investigated both in steady state and transient scenarios. The investigated transients consisted in partial assembly flow blockage. It represents both a severe accident initiator in SFRs, as seen in the partial fuel meltdown in the

Unit 1 of the E. Fermi Atomic Power Plant (Michigan, U.S.A.) (Bertini, 1980) and a suitable case to highlight the impact of 3-D effects due to the large radial thermal gradients and sodium flow through windows that arise. The analysis is inclusive of both a mesh convergence analysis and a verification effort of the newly developed models for treating the inter-assembly heat transfer. The paper is organized as follows In Section 2, the theoretical foundations of a porous medium based coarsemesh approach are presented alongside the development of a novel model for heat transfer across the assembly-wrapper for the GeN-Foam code. In Section 3, the modelling of the assembly bundle and the blockage investigation procedure are presented. In Section 4, the results of the mesh convergence study and verification against analytical results are presented. In Section 5, the results of the actual parametric analysis of the thermal-hydraulic impact of the features under investigation are presented. In Section 6, conclusion and future work directions are drawn. 2. Methodology and further developments 2.1. A porous medium based coarse-mesh methodology A porous medium based coarse-mesh approach and a sub-channel approach are related as both rely on correlations for the closure of momentum and heat transfer terms in the governing equations, providing comparable degrees of accuracy. However, the former can be seen as a generalization of the latter. In sub-channel codes, the bundle geometry needs to be explicitly resolved, and the modelling capabilities of various sub-channel codes are typically hard-coded for specific assembly geometries. In contrast, by borrowing the mathematical tools commonly employed in porous media theory (Vafai, 2005), namely volume averaging techniques, a coarse-mesh methodology enables the description of any structure of interest with general unstructured meshes and without the need for an explicit geometric representation of the so-called sub-scale structures. By sub-scale one means any structure whose characteristic dimension is smaller than the characteristic dimension of the mesh cells employed to describe the system. This loss of detail is compensated with correlations for momentum exchange (i.e. pressure drops) and heat transfer (which depend on the bundle geometry and should be chosen accordingly). These correlations are not forceful additions to a coarse-mesh description of a system of interest. Instead, these naturally arise as additional terms in the volume averaged RANS equation that happen to require adequate closure. It should be reminded that the volume averaging process comes at the cost of a possible loss of accuracy. From a mathematical perspective, a volume averaging technique provides a macroscopic system description starting from a set of locally governing equations. By “locally” one means at the sub-scale level (e.g. coolant in between fuel pins in a reactor assembly), while by macroscopic one means at a coarse-mesh level. For clarity, an example is provided in Fig. 2.1. Let us consider a volume V which encompasses both a fluid phase

Fig. 2.1. Left: geometry of a fast reactor assembly and the fuel pin lattice. Right: a possible coarse-mesh description of the assembly. Each cell encompasses both the fluid and the sub-scale structure (i.e. fuel pins in this case). 2

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

(denoted by ) and a solid phase (denoted by , e.g. the sub-scale structure), and any scalar or vector field * locally defined in the fluid phase. The intrinsic and superficial average operators are defined as, respectively: i

1 V

=

=

1 V

dV

V

p

(2.2)

The star symbol indicates that the quantity is local. The intrinsic average represents the physical average of the selected quantity over the volume occupied by the fluid within V, while the superficial average represents its average over the entire volume V. This distinction comes into play as intensive physical quantities (with the notable exception of velocity) are typically described via their intrinsic average within a coarse-mesh cell, while extensive physical quantities are described in terms of their superficial average (Vafai, 2005). A quantity that relates the two averages is the porosity . It is defined as fluid volume over total cell volume V. It can be noticed that:

=

t

+

t

·( u ) = 0

( u)+

(2.4)

·( u

u)=

p +

·S +

g

(2.5)

For generality, the form of the deviatoric stress tensor S is not specified. Due to the additivity of the averaging operators, applying the superficial average operator to the equations consists in superficially averaging each term of the equations. To this regard, it can be show that the volume average of temporal or spatial derivative operators results in:

=

t =

·

=

1 V

t +

·

1 V

+

1 V

V,

V,

V,

us · ns dS ns dS

·ns dS

=

p

+

· S

+

1 V

V,

( p I + S)· ns dS

(2.9)

( u) +

·

1

u

u =

( p) +

·(µt u) + F , +

g

(2.10)

where the variables represent the intrinsic volume averages of the corresponding local variables over each mesh cell, with the exception of the velocity. The velocity field in the equations is the superficial average of the corresponding local velocity field in each mesh cell, also known as the Darcy velocity field. That is the velocity that the fluid would have if it occupied the entire volume of the computational mesh unit. The density is maintained in the momentum equation as an incompressible flow does not necessarily imply an incompressible fluid, as density changes might still arise in the domain (e.g. due to temperature gradients). These equations closely resemble their microscopic counterpart except for the presence of the additional momentum source/sink term F , on the right hand side of the momentum equation, and the presence of the porosity term . What should be noticed is that if the porosity approaches unity (i.e. towards the limit of no subscale structures in a mesh cell, also called a “clear fluid” condition), these equations revert back their standard RANS formulation employed in traditional CFD. From a computational perspective, this implies that the methodology can be used to represent an entire computational domain that is inclusive of both porous coarse-mesh regions and clear fine-mesh regions without the need to segregate the matrix description of the domain. For this reason, this approach is also referred to as unified fine-coarse-mesh approach, yet for simplicity it will be referred to as coarse-mesh approach for the rest of the present document. In terms of the energy equation, while the volume averaging of the microscopic energy balance gives rise to an equation that is analogous to Eq. (2.10), the treatment of the thermal diffusivity would require some attention. As a matter of fact, the topology of the porous medium will alter the effective thermal diffusivity of the fluid, and this is typically quantified by what is commonly referred to as tortuosity of the porous medium. In principle, the tortuosity is a rank 2 tensor quantifying the deformation of the path along which a quantity diffuses (e.g. internal energy in this case), and is meant to adjust the effective thermal diffusivity in the energy equation.

In order to obtain the macroscopic governing equations, the local RANS equations are superficially averaged over each mesh cell. For simplicity, only the continuity and the momentum equations are considered. Their local formulation is:

t

·S

·u = 0

(2.3)

i

+

The additional term thus consists in the surface integral of the total stress tensor at the interface, which physically represents the momentum exchange per unit time per unit volume between the two phases, and will be referred to as F , . From the perspective of an additional energy equation, it can be shown that the volume average of the heat diffusive term will give rise to not necessarily null contributions from the surface integrals at the phase interfaces. This in turn represents the energy exchange between the two phases, and requires closure as well. The volume averaged RANS equations, can thus be obtained at a coarse-mesh level. For an incompressible flow, these consist in:

(2.1)

dV

V

rise to generally non-null surface integrals. In particular:

(2.6) (2.7) (2.8)

The integrals span over the interface between the two phases V , , and the interface normals ns are oriented from the solid phase outwards to the fluid phase . The volume average of a time derivative is presented in its most general form for a non-stationary interface moving at velocity us , yet this term is zero for any case of interest in nuclear applications. Without recurring to the effective yet lengthy derivation offered by Whitaker (1985), the presence of surface integrals is qualitatively justified. A volume average of a divergence or gradient consists in a volume integral of these quantities, which can be transformed to a surface integral via the Gauss-Green theorem. However, within the averaging volume (i.e. the coarse-mesh cell from a computational perspective) the integrated fluid field has discontinuities at the phase boundaries, so that the surface integral will need to account for the phase boundaries internal to the averaging volume. It is these integrals that represent the additional terms that require closure. As a matter of fact, let us consider the volume averaged momentum equation. The volume average of the advective term will have a null contribution from the surface integrals at the phase interfaces due to the local no-slip condition at the phase interface (i.e. fluid velocity is null at the interfaces with solid boundaries). Conversely, the volume average of the pressure gradient and the divergence of the deviatoric stress tensor give

2.2. Implementation in GeN-Foam The goal of the following sub-sections is to illustrate some features and implementation choices of the coarse-mesh methodology in GeNFoam. For a thorough discussion on the matter, refer to Fiorina et al. (2015). As it was previously discussed, the volume average of the diffusive terms in a governing equation will give rise to additional surface integrals to be evaluated at the boundaries between the phases within the averaging volume. This represents an inter-phase transfer term, either for momentum or energy and should be modelled adequately. 3

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

2.2.1. Momentum transfer term The inter-phase (i.e. drag between fluid and sub-scale) momentum transfer term is proportional to fluid velocity in a manner that is usually in between linear and quadratic proportionality. Thus, GeN-Foam represents such term (i.e. F , in Eq. (2.10)) in the following manner:

computed via correlations for the turbulent intensity It and turbulent length scale Lt for channel flows (ANSYS FLUENT, 2019):

keq =

(2.11)

F , = D (u)· u

3

eq

where D (u) is a second rank drag tensor that captures any higher order dependence of drag on velocity other than the linear one. This allows to include the drag term in a computationally semi-implicit manner in the momentum equation, favouring stability. For choices of a local frame of reference that are aligned with the main flow direction, the drag tensor can be simplified to a diagonal tensor whose components can be replaced with established correlations for pressure drop computation depending on the sub-scale geometry. Because of this, GeN-Foam offers the option to specify a local reference frame for any region of interest in the computational domain. On top of that, the OpenFOAM framework of the code allows for simplified implementation of novel pressure drop correlations.

h=

Tf )

Nu = ANu + BNu PeCNu

(2.15)

3

k2 Lt

(2.16) 0.125

Lt = 0.07 Dh

(2.17) (2.18)

The present turbulent approach in the porous medium thus relies on standard experimental correlations for the estimation of the turbulent parameters. For simplicity, and due to the lack of detailed experimental data, standard correlations for channel flow have been considered as suitable for the scope of this work. However, the effects of wire-induced mixing on the thermal and momentum diffusivity could in principle be captured by fine tuning the turbulence parameters in order to fit experimental data (as done e.g. by Khan et al. (1975) in past porous medium approach investigations). 2.3. Boundary conditions for coarse mesh models 2.3.1. Wrapper heat transfer modelling In SFRs, the fuel assemblies are encased in so-called wrappers and are separated by a sodium gap a few millimeters in width. These can be modelled via internal boundaries (called baffles) that physically separate the flow in assemblies and gap, and that take into account the thermal resistance of the wrapper. However, from a coarse-mesh point of view, the unit mesh cell size will typically be significantly larger than the size of the boundary layer. An adequate wall model is then needed for computing the additional heat resistance of the latter. To understand the equations that are being used to simulate heat transfer through wrappers in GeN-Foam, let us consider a baffle with a physical thickness b , a thermal conductivity kb and let us denote the two sides of the baffle with + and − superscripts. The temperatures of two opposing baffle parcels will be denoted with Tb+ and Tb , while the cell center temperatures of the cells adjacent to these baffle parcels will be denoted with Tc+ and Tc . A representation of this is provided in Fig. 2.2. However, it should be stressed that for ease of geometry and mesh generation, the baffle thickness is not explicitly geometrically resolved in the computational model. For this reason, the two sides of the baffle have the same surface area. By neglecting the thermal inertia of the

(2.12)

Nu kf Dh

= Cµ4

It = 0.16 Re

2.2.2. Heat source term A sub-scale model can be prescribed in any porous region and associated to a temperature Tss and a heat capacity. When the fluid equations are solved for, a volumetric heat source (i.e. the one that arises from coarse-mesh modelling that requires closure) is computed as:

Q = Ass h (Tss

3 (|u|It )2 2

(2.13) (2.14)

where Ass is the average surface area of the subscale structure per unit volume, Tss and Tf are the subscale and fluid temperatures (within the mesh cell) respectively. The heat transfer coefficient h is computed via adequate correlations for the Nusselt number in their typical form for liquid metals in pin bundles (Eq. (2.14)). Different set of coefficients ANu , BNu and CNu can be specified for the laminar or turbulent flow regimes. Given the nuclear vocation of the code, a sub-scale fuel model can be prescribed in any porous region that represents fuel. In particular, within each mesh cell, the model consists in a finite difference 1D radial heat diffusion problem through fuel, gap and cladding. A volumetric power source in the fuel can be either prescribed for singlephysics calculations or provided by the neutronics solver. When the 1-D heat equation is being solved for in the fuel, a convective boundary condition is applied to the cladding with the available fluid temperature and heat transfer coefficient values. The solution is thus obtained via operator-splitting and an implicit coupling between the fluid-dynamics and fuel model is achieved via fixed-point iterations. Both the 1-D meshing and geometric properties of the fuel pins are user selectable. No axial heat conduction is modelled due to the typically smaller axial temperature gradients in oxide fuels. 2.2.3. Turbulence modelling While the prediction of turbulent fields in porous media enjoys a level of interest (Chandesris et al., 2006; Nakayama and Kuwahara, 1999), it remains of limited engineering applicability. However, it should be recalled that the goal of a turbulence model is to adequately close the RANS equations so that, from an engineering perspective, effective fluid viscosity and thermal diffusivity are can be captured satisfactorily. To this regard, While GeN-Foam relies on a traditional k (Launder and Sharma, 1974) approach for clear fluid regions, the k and equations are not solved for in coarse-mesh regions. Rather, the values of k and are forced to exponentially converge, along the flow length, to equilibrium values. These equilibrium values keq, eq are

Fig. 2.2. Schematic of two computational cells at the opposing sides of a thermally conductive baffle. Cell centers are marked with full circles while face centers at cell boundaries are marked with half circles. The heat transfer coefficients between relevant positions in the geometry are reported. From a geometric modelling perspective, the baffle has no thickness but the heat transfer is treated as if it had a thickness b . 4

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

baffle (a reasonable assumption if the absolute heat capacity of the wrapper is small compared to that of the sodium in the domain), the baffle surface temperatures can be computed by enforcing heat flux conservation, namely that the heat per unit area that flows into one side of the baffle is instantaneously the same as the heat per unit area flowing out of the baffle. Thus:

1 b

kb

+

1 h+

1 b

kb

+

1 h

(Tc+

Tb ) = h (Tb

Tc )

(Tc

Tb+) = h+ (Tb+

Tc+)

OpenFOAM and applied to the windows. Currently it consists of an explicit correction to the pressure term, so that a number of PIMPLE iterations are required for the momentum and pressure equations to reach convergence. Stability was not found to be an issue due to the relatively low velocity magnitudes that characterize the system. 3. Modelling and investigation procedure

(2.19)

3.1. System modelling Three separate geometric models are investigated, based on the data of the ESFR-SMART core assemblies (Appendix A). These models are: 1) a single assembly with thermally adiabatic wrappers; 2) a bundle of seven assemblies with thermally conductive, window-less wrappers, separated by an inter-assembly gap; 3) a bundle of seven assemblies with thermally conductive wrappers with windows, separated by an inter-assembly gap. These models will be referred to as “adiabatic”, “windowless” and “windowed” respectively from now on. The adiabatic and windowless models span over the axial length of the fuel region within the assemblies. The windowed model extends up to the assembly outlet. This modelling extension is necessary as sodium flow through windows is governed by the pressure gradients between the assemblies and the gap. This in turn depends on the entire axial structure of the assembly, as a common pressure is only shared at the outlet. An overview of the geometry of the models is presented in Figs. 3.1 and 3.2. The model is porous medium based, meaning that the fuel pins are not geometrically resolved, and the assembly properties are volume averaged. In particular, the porous properties of the assemblies were obtained by averaging the pins over the inner assembly volume (i.e. exclusive of the wrapper), while the gap properties were obtained by averaging the wrapper thickness over the volume of the wrappers and the inter-assembly gap volume. Relevant geometric data, porous medium modelling data, friction factor and Nusselt number correlations with the related parameters, if applicable, are reported in Table 3.1. In particular, the data for each of the axial regions of the windowed model is reported, while the adiabatic and windowless models share the data of the fuel axial region of the windowed model. With reference to Table 3.1, a few remarks should be made. The gap height was not reported as it corresponds to the total model height, which varies between the three models. Two values of the hydraulic diameter for the fuel region are reported. The “bulk” value is employed for the computation of the friction factors and pin-sodium heat transfer coefficients (via Eq. (2.15)), while the “wall” value is used for sodiumwrapper heat transfer coefficient calculations. With regard to friction factor correlations, a Blasius correlation is used to model friction in the intra-assembly plena, heads and in the inter-assembly gap. For the pinoccupied regions, a Rehme correlation (Rehme, 1973) with adequate geometric parameters was employed for the calculation of the momentum sink parallel to the bundle direction, while a Gunter-Shaw correlation (Gunter and Shaw, 1945) was employed for the momentum sink due to cross-flow in the direction perpendicular to the pins. Two sets of parameters for the Nusselt calculations are employed throughout the model, one for the laminar regime (i.e. for Re < Relam trans ) and one for the turbulent regime (i.e. for Re > Retrans turb ). The coefficients were obtained from the proposed correlation by Mikityuk (Mikityuk, 2009) with respect to the present fuel bundle geometry for the turbulent regime. For the laminar regime, the Nusselt number is constant. For the transition regime, the Nusselt number is obtained via linear interpolation between the two Nusselt numbers computed at the transition regime boundaries. While the use the same correlation (devised for pin bundles) for the Nusselt number for both the pin-fluid heat transfer and the wrapper wall-fluid heat transfer is debatable, it should be noted that this has little impact on the results of the present work. As a matter of fact, the overall heat transfer coefficient between the intraassembly sodium and the sodium in the assembly gap will be limited by the thermal resistance of the wrapper wall (i.e. 4.5 mm of 20 W / m /K

(2.20)

The equations represent an instantaneous heat flux balance, given that: 1) the heat transfer between a cell and the adjacent baffle parcel on the + side is characterised by the heat transfer coefficient h+; 2) the heat transfer between the two opposing baffle parcels on the + and − side respectively by the heat transfer coefficient b ; 3) the heat transfer kb

between the baffle parcel and the adjacent cell on the − side is characterised by heat transfer coefficient h . These heat transfer coefficients are computed according to Eqs. (2.13) and (2.14). The coefficients ANu , BNu and CNu user selectable and different sets of coefficients can be specified for the laminar and turbulent flow regimes, with linear interpolation in the transition regime. Within each time step, the procedure is to: 1) solve for the energy equation in the domain provided guess values or previous iteration values for the baffle temperatures Tb+, Tb as fixed values; 2) update baffle temperatures based on Eqs. (2.19) and (2.20); 3) repeat starting from step 1) until convergence. More specifically, these iterations are performed at the level of the outer iteration loop of the PIMPLE algorithm (Holzmann, 2018), employed by the GeN-Foam code. In principle such boundary could be rendered implicit by directly manipulating the discretised energy equation matrix coefficients corresponding to mesh cells at the opposing sides of the wrapper (thus sparing the need for the iterative process within each time step). However, the overhead in terms of energy equation iterations that the boundary adds (the total iteration count being the range from 6 to 10) is not a concern as of now. Developments in this direction are nonetheless planned. It should be noted that the present heat transfer model only treats heat transfer perpendicular to the wrapper surface. 2.3.2. Pressure baffle Consider two adjacent regions of the computational domain characterized by some macroscopic averaged parameters but whose subscale structures are geometrically different. Given a fluid flow crossing both regions, an additional irreversible pressure jump term is physically expected at the interface because of the actual change of the sub scale structure geometry (causing a non-isoentropic change in fluid state). For the practical case under investigation of an ESFR assembly, this scenario is identified at the wrapper windows. In particular, the pressure jump is modelled in a similar manner as for orifices and pipe bends, that is:

p=

1 |u|2 2

(2.21)

where is the head loss coefficient. As a limiting case (i.e. maximum pressure drop), the head loss coefficient for the present application was assumed to be equal to that associated with a double right-angled mitered pipe bend. This would be the pressure loss if the fluid e.g. entering the assembly from the gap were forced to bend at a right angle with respect to the main flow direction. Additional irreversible pressure losses are expected within the assembly itself at the interface between regions characterized by sharp changes in the porosity (e.g. the pin bundle and the intra-assembly plenum). However, these losses were found to be negligible in magnitude. Such a boundary condition for the non-hydrostatic component of the pressure was developed in 5

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

Fig. 3.2. Detail of the assembly wrapper: a) without windows; b) with windows. Please notice that the wrappers in the windowed model extend up to the assembly outlet. However, for a clearer view of the windows, the wrappers have been cut at the end of the fuel region.

hydraulic diameter. A head loss coefficient of = 1.17 was assumed for the pressure drop across the windows. In particular, this value was obtained by averaging the Reynolds-dependent correlation for head loss coefficients of right angle miter bends devised by Al-Tameemi and Ricco (1975) over the range of Reynolds numbers of the investigated cases. The magnitude of the pressure drop is then computed according to 2.21, with appropriate directionality taken into account and the velocity being the real fluid velocity and not the Darcy velocity. For the simulation of the different investigated scenarios, fixed value boundary conditions are imposed at the inlets of central assembly, neighbouring assemblies and gap for the velocity and a common pressure is set at all the outlets. An axially cosine shaped, radially uniform fuel power profile is set across the fuel region in each assembly to model internal power production.

Fig. 3.1. Side view and a quarter of top view for: a) the single assembly model (i.e adiabatic); b) the seven assembly bundle with thermally conductive wrappers without windows (i.e windowless); c) the seven assembly bundle with thermally conductive wrappers with windows (i.e windowed); please notice that the windows cannot be seen from the side view as they are shrouded by the half-gap around the assemblies. Nonetheless, the location of the windows behind the half-gap mesh is highlighted in yellow. The windowed model consists of five axial regions: I) fuel; II) plugs; III) intra-assembly plenum; IV) upper shield; V) assembly head. d) Detail of the inter-assembly gap meshing. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

steel). In particular, even the worst case estimates for the fluid-wrapper 3) result in a heat transfer coefficient that is Nusselt (e.g. for a Nu one order of magnitude larger than the heat transfer coefficient of the wrapper itself. With regards to the choice for the tortuosity of the porous medium (see Section 2.1), this was chosen to be equal to the volumetric porosity of the pin bundle. While other investigations suggest the use of the minimum lateral surface porosity of the bundle (e.g. (Khan et al., 1975; Yu et al., 2015)), this option was deemed too conservative in our scenario, as said porosity is considerably smaller than the volumetric porosity due to the small pitch-diameter ratio of the bundle. Fine-tuning the turbulence model to fit currently unavailable experimental data was out of scope for the current work. With regard to wrapper windows, these are opening with a hydraulic diameter of 60 mm whose center lies at 168 mm below the top of the fuel. It is recalled that these opening are located close to the top of the fuel so to prevent, in boiling scenarios, downward axial vapour propagation towards core regions characterized by a locally positive void reactivity coefficient. While the windows are circular in shape, for modelling simplicity these were represented as square openings with the same

3.2. Investigation methodology Each simulated case is defined by four parameters, namely: 1) reactor power expressed as a fraction fP of the power at steady state operation; 2) total sodium velocity in the neighbouring assemblies expressed as a fraction fN of the reference sodium velocity at steady state operation; 3) sodium velocity in the central assembly expressed as a fraction fC of the sodium velocity in the surrounding unblocked assemblies; 4) sodium velocity in the inter-assembly gap expressed as a fraction fG of the sodium velocity in the surrounding unblocked assemblies. A parametric approach with respect to this quantity is justified by the lack of experimental data regarding inter-assembly gap velocity magnitudes. A list of values of interest for each of these quantities is made. For every combination of these values, for each of the three models (i.e. adiabatic, windowless, windowed), simulations are performed in the following way. First, a simulation with no blockage is performed until (and if) a steady-state is reached, given the other parameters for reactor power, unblocked assembly sodium velocity and gap sodium velocity fractions (the gap is not modelled in the 6

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

Table 3.1 Geometric and porous medium modelling data for the axial regions of the windowed model. The adiabatic and windowless models share the same data of the fuel region of the windowed model. Quantity

Height [mm] Width [mm] [] Dh [mm]

Ac

m2 m3

Friction factor correlation

Nusselt correlation

Re lam Re trans

trans turb

Value in region Fuel

Plugs

Plenum

Shield

Head

Gap

864.29 196.35 0.268 3.35 (bulk) 4.36 (edge) 272.08

68.00 196.35 0.268 3.35

600.00 196.35 1.000 113.36

576.00 196.35 0.280 3.37

230.00 196.35 1.000 113.36

– 13.50 0.362 9.10

0

0

0

0

0

Rehme & Gunter-Shaw 271 N [] 10.67 Dp [mm]

Rehme & Gunter-Shaw 271 N [] 10.67 Dp [mm]

Blasius

Rehme & Gunter-Shaw 19 N [] 40.12 Dp [mm]

Blasius

Blasius

Pp [mm]

Pp [mm]

D w [mm] Hw [mm] Value Laminar regime ANu BNu CNu Turbulent regime ANu BNu CNu 1000 2300 1.17

11.67

1.00 225.00

D w [mm] Hw [mm]

11.67

Pp [mm]

1.00 225.00

D w [mm] Hw [mm]

43.15

1.00 225.00

3.521 0 0 3.521 0.014 0.770

adiabatic case for obvious reasons). Afterwards, an instantaneous blockage is modelled in the form of a sudden sodium velocity drop in the central assembly, and the simulation is performed until a new steady state or if boiling temperature are reached. If a new steady state is reached, the outlet bulk temperatures (i.e integral outlet temperature weighted by the enthalpy outflow per unit time) of the gap and central assembly are estimated. If boiling temperature is reached, only the time it took to reach boiling is registered and simulations are terminated. With regard to this work, the investigated parameters consist in: 1) 100%, 5% reactor power; 2) 100%, 5% sodium velocity in surrounding assemblies; 3) sodium velocity in central assembly from 0% (i.e total blockage) to 60% in steps of 10% of the neighbouring assemblies sodium velocity; 4) sodium velocity in gap from 0% to 100% in steps of 20% of the neighbouring assemblies sodium velocity. In particular, the investigated power levels were selected so that results characteristic of varying degrees of blockages at both full power ( fP = 100%) and at power levels compatible with core status right after shutdown ( fP = 5%, e.g. a protected transient following blockage detection) could be obtained. Nominal assembly parameters correspond to a power of 7.143 MW and a flow rate of 44.7 l/s. This corresponds to a sodium velocity of 5.00 m/s through the active core pin region.

TG (x ) x TA (x ) x

= =

S h w L

|uG | SG cp

(TA (x )

S h w L

|uA | SA cp

(TA (x )

TG (x )) TG (x ))

(4.1)

Provided with the thermo-physical properties of the fluid (density and heat capacity cp ), channel geometric data (shared surface area Sw and total channel length L) and inlet boundary conditions for the velocity uA,0, u G,0 and temperature TA,0, TG,0 , a general solution can be obtained in the form:

TG (x ) = c1 v1,1 e

1x

+ c2 v2,1 e

2x

TA (x ) = c1 v1,2 e

1x

+ c2 v2,2 e

2x

(4.2)

were i is the i th eigenvalue of the matrix representation of system 4.1, vi, j is the j th component of the eigenvector associated with the i th eigenvalue, the coefficients ci are to be obtained via the imposition of boundary conditions. The sodium properties and geometric data used for the verification are reported in Table 4.1. The GeN-Foam model of the channel alongside the obtained results and discrepancies with respect to the analytic model are presented in Fig. 4.1. To be consistent with the assumptions made to derive the analytic solution, only one radial mesh cell was used both in the assembly and the gap, while a total of 75 cells were employed for the axial discretisation. The verification is deemed acceptable given the magnitude of the relative discrepancies between the two axial temperature profiles in the gap and in the assembly.

4. Verification and mesh convergence study 4.1. Heat transfer model verification The newly developed heat transfer model was verified against an analytical solution prior to the actual blockage investigation. The analytic case is a steady state 1-D problem consisting of parallel fluid flow in two channels with fixed inlet temperatures that share a certain surface area through which heat can be exchanged. Let us denote the channels with the labels A and G and let us assume that these are parallel to a certain axis x . If there is no internal power production within the channels, the spatial temperature distribution is governed by:

4.2. Mesh convergence study 4.2.1. Theoretical remarks A thorough mesh convergence analysis was motivated by the specific challenges that coarse meshes pose. While a coarse-mesh has less room for improvement in terms of actual mesh-independence of the results compared to standard fine-mesh approaches, it can still be a valid choice provided that the obtained results are in the asymptotic 7

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

0.72

solution. To avoid large discretization error one should be able to make an estimate on the error bounds obtained on the quantities of interest to be computed by the simulation. For this estimate to be possible, however, the results need to be in what is referred to as asymptotic convergence range. This is typically not assured for results obtained in a coarse-mesh domain (Roache, 1994). For a result of interest to be in the asymptotic convergence range, it needs to monotonically approach a definite value as the mesh is further refined. This behaviour enables a meaningful extrapolation of the results at an ideally infinitesimal mesh size, such as generalized Richardson extrapolation (Sidi, 1997), which in turn allows for error bounds to be estimated. Let us consider a certain output quantity of interest f. If fe is the extrapolated value at an ideally null mesh spacing and fh is the value obtained at a certain mesh spacing h, the error Eh can be estimated as:

820.00

Eh = fh

1.34

In particular, the obtained result fh is in the asymptotic convergence range if the error Eh can be expressed as a polynomial of the grid spacing whose lowest order (i.e. of order p) coefficient C doest not depend on h. The big O notation is used to represent higher order terms that become negligible as the mesh spacing approaches infinitesimal values. If the results are in the asymptotic convergence range, given three output values obtained respectively on progressively

Table 4.1 Geometric data, sodium properties and boundary conditions used for the verification. Quantity

Value

L [cm]

75.00 5101.34

Sw [cm2] h

W m2 K

kg m3

cp

J kg K

TG, 0 [K]

SG [cm2] |u G |

m s

TA, 0 [K]

SA [cm2]

|uA |

m s

3693.72 841.06 1268.53 668.00 47.49

338.83

convergence range (Roache, 1994). This enables the extrapolation of the results at null grid spacing via techniques such as Generalized Richardson Extrapolation (Sidi, 1997) and thus enables the estimation of expected error bands. When selecting a mesh for a case of interest in a coarse-mesh context, it should be recalled that: 1) an overly fine mesh would defy the purpose of a coarse-mesh approach; 2) an overly coarse mesh might lead to unacceptable discretization errors in the obtained

(4.3)

fe = C h p + O (p > p)

h

h

refined meshes by a factor r = h i = hi + 1 , the convergence order p and i+1 i+2 the extrapolated value f can be obtained as (via Richardson extrapolation):

Fig. 4.1. Geometric model of the channel in GeN-Foam (left) and results of the verification (right). The black curve represents the temperature profile obtained at steady state by GeN-Foam, while the red curve represents the relative discrepancy between the solutions obtained by GeN-Foam and the analytic ones.

8

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

by a coarse mesh fi is in the asymptotic convergence range if:

GCIi

GCIi + 1 r p

(4.10)

This can be computed as long as the output parameters fi + 1 , fi + 2 obtained at subsequently refined meshes by a factor r are available. Thus, to assess whether results are in the asymptotic convergence range, the following ratio is evaluated and checked whether it is within unity:

ACR =

p=

(

fi

fi + 1

fi + 1

fi + 2

)

fe

fi + (fi + 1

4.2.2. Results Two cases were investigated, namely: 1) fP = 5%, fN = fP = 100%, fN = 100%, fC = 50%, fG = 100%; 2) 5%, fC = 50%, fG = 100% . These two cases will be referred to as “full power” and “low power”, respectively. The quantity for which convergence properties were assessed was the outlet bulk temperature of the central assembly and the gap. In particular, with regard to the results concerning the GCIs, subscripts 1 and 2 will be used to denote the GCIs related to the coarse and fine meshes respectively. Windowless model. For both the full and low-power cases, two different mesh topologies were investigated, reported in Fig. 4.2. Referring to said figure, these meshes will be referred to as triangular mesh and web mesh. The reason for investigating these topologies is that while the triangular mesh is among the most popular choices for the investigation of fast reactor assemblies, it might not be the most suitable option due to the added numerical diffusivity in this scenario. This comes from the fact that many cell faces will not be orthogonal to thermal gradients. The web mesh has been investigated to try to address this issue. The radial divisions within an assembly are such that each cell has the same overall volume. For each mesh topology, different meshes obtained by varying the number of nodes in the characteristic directions were investigated: 1) nA number of nodes in the axial direction; 2) nR number of nodes in the assembly radial direction; 3) nS number of nodes per assembly side (only for web meshes); 4) nG number of nodes in the half-gap radial direction. Any mesh for the windowless model can be described in terms of these four parameters. The mesh that was found to be suitable for the scope of this work is a web mesh characterized by nA = 20, nR = 4, nS = 1, nG = 4 . Convergence properties were investigated separately for radial and axial mesh refinements. It is recalled that in order to assess the convergence properties, the results need to be obtained on two further refined meshes. For the case of radial refinement, this consisted in subsequently doubling both the number or radial nodes in the assembly nR and half-gap nG , while maintaining the same number of nodes in other directions. For the case of axial refinement, this consisted in subsequently doubling nA while maintaining the same number of nodes in the other directions. The mesh convergence properties of the outlet bulk temperatures for the central assembly TC and gap TG for the low power case are provided in Table 4.2. The results are found to be in the asymptotic convergence range 1) for both the triangular and web meshes and all the re( ACR finement types. Convergence properties related to central assembly temperature TC for axial mesh refinement are not reported as these were found to be already at geometric convergence (i.e. < 0.1 K difference between the coarse and the extrapolated values). For a given set of radial mesh parameters, further axial refinements do not have a significant impact on the solution. This is expected, as the thermal

(4.4)

ln (r ) fi )

rp rp

(4.5)

1

An error estimator on the coarse solution fi can thus be obtained. Various estimators for the error on the coarse mesh solution can be defined: i

fi

=

fi + 1 (4.6)

fi fi

ei =

fe

=

fi

i

rp rp

(4.7)

1

Neither of these should be interpreted as error bounds on the coarse solution fi , rather as error bands, i.e. a tolerance on the accuracy of the solutions, that may as well be exceeded. In the experience of Roache (1994), the most commonly reported error estimator in mesh convergence studies in the one represented by Eq. (4.6). However, this estimator contains no information of the actual order of convergence of the solution on its own, which is necessary to assess if the solution is in the asymptotic convergence range. This information is instead present in the relative error estimator ei . To address this, Roache proposes the use of a so-called Grid Convergence Index (GCI) defined so that: i

rp rp

1

= GCIi

p r0 0 p0 r0

(4.8)

1

The GCI thus represents what would be the equivalent error estimator in the form of Eq. (4.6) if the results were obtained by a p0 -order accurate method for a mesh refinement ratio r0 . Typical values suggested by Roache for the comparison are r0 = 2 and p0 = 2 , thus consisting in a comparison against the uncertainties associated with the convergence of a second order accurate method with mesh halving. Thus, from Eqs. (4.6)–(4.8) it follows that:

GCIi =

3 3 ei = 4 4

fi

fi + 1 fi

rp rp

1

(4.11)

To summarize, the mesh convergence effort is aimed at establishing whether or not results are in the asymptotic convergence range via the estimation of the ACR and the subsequent estimation of representative expected error bands on the obtained results. This process is furthermore meant to establish which mesh topology is best suited for the spatial discretization of the domain given the nature of the investigated cases.

Fig. 4.2. Examples of half cross-sections of the investigated mesh topologies: (left) triangular mesh; (right) web mesh. The key quantities that characterise each mesh are reported. These consists in the number of computational nodes in each reported direction. Each mesh is further characterized by the number of nodes in the axial direction nA , not reported in the figure.

ln

GCIi GCIi + 1 r p

(4.9)

It should be noted that these formulations differ by those of Roache by a factor 1p0 as this is done from the perspective of error estimation r0

for the coarse solution fi and not for the fine solution fi + 1, yet the approach is the same. It can be shown that the output of interest obtained 9

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

Table 4.2 Mesh convergence results for the outlet bulk temperature of the central assembly TC and gap TG for the low power case. The results are presented for radial and axial refinements starting from the results obtained by a web mesh and a triangular mesh with the same parameters: nR = 4, nS = 1, nG = 4, nA = 20 . The quantity p represents the order of convergence as defined in Eq. 4.4. The convergence results related to the axial refinement of the mesh on TC are not reported because the obtained values for TC were found to be already at geometric convergence (i.e. <0.1 K difference between the coarse and the extrapolated values) for both the web and triangular mesh.

TC

TG

Mesh

Value [K]

Refinement

p[]

GCI1 [%]

GCI2 [%]

ACR [%]

Ex. val. [K]

Err. bnd. [K]

Tri

894.17

1.27

0.524

0.218

6.30

894.87

1.32

0.429

0.172

100.170 – 100.137 –

900.44

Web

Radial Axial Radial Axial

900.01

5.15

Tri

771.64

Web

771.10

Radial Axial Radial Axial

0.70 0.67 0.71 0.56

1.000 0.054 0.817 0.051

0.614 0.034 0.500 0.035

99.685 100.017 99.741 100.015

761.44 772.20 762.74 771.63

10.12 0.56 8.31 0.5

gradients in the axial direction are much smaller than those in the radial direction. Furthermore, the web mesh is characterized by better mesh convergence properties in terms of the estimated error bands on the obtained results. This was also expected due to the previous remarks on the reduced numerical diffusivity that the web mesh provides. The same trends in the mesh convergence analysis were found in the full power case. In particular, the results were still consistently in the asymptotic convergence range and the web mesh proved to be a more desirable choice with respect to the triangular mesh. The associated uncertainty bands were approximately twice as smaller when compared to the low power scenario. For consistency, the same mesh was employed for the adiabatic case. Windowed model. The mesh choice for the windowed model is defined by a larger number of parameters with respect to the windowless model as there are multiple axial regions (refer to Fig. 3.1c). The mesh convergence properties of a windowed mesh with similar mesh parameters as those selected for the windowless model were assessed. In particular, the investigated windowed mesh consisted in a web mesh with nR = 4, nG = 4 (i.e. the same as for the windowless model), nS = 4 (the bare minimum required to model windows of the correct size) and a variable number of mesh nodes depending on the axial zone. Separate axial and radial refinement studies were performed on the mesh for the low power and high power cases. Results for the web mesh only (the triangular mesh was not investigated) for the high power case are presented in Table 4.3. All of the results are in the asymptotic convergence range, and the same was observed for the low power case. The associated uncertainty bands are of the same order of magnitude as for the windowless case. However, the relative importance of axial refinement with respect to radial refinement is increased in this case, as testified by the reported uncertainty bands. The same was found for the low power case. This was expected, as the outlet bulk temperatures are no longer shaped exclusively by the heat transfer process across the wrapper, but also by the 3-D mixing processes that occur at the windows. As an example, the differences in flow resolution around wrapper windows between the coarsest and the most refined (both axially and radially) case is presented in Fig. 4.3. In conclusion, the selected web mesh that is radially

identical to and axially comparable to the one used for the windowless model was found to be suitable for the windowed model as well. 5. Assessment of inter-assembly gap and wrapper window behaviour The investigation was performed parametrically for two scenarios: 1) full power and flow conditions corresponding to nominal operation values ( fP = 100%, fN = 100%), which will be referred to as “full power” case; 2) low power and flow conditions compatible with core status in natural circulation regime after a protected transient ( fP = 5%, fN = 5%), which will be referred to as “low power” case. For each of the two scenarios, simulations were performed for different combinations of the central assembly sodium velocity (to model blockage) and inter-assembly gap sodium velocity. Each simulation is performed on all three models, “adiabatic”, “windowless” and “windowed”. To better visualize the thermal-hydraulics of the system, results obtained by the simulation of one of the multiple investigated cases is reported in Fig. 5.1. 5.1. Full power, full flow A comparison of the central assembly outlet bulk temperatures, for varying degrees of central assembly sodium velocity and gap velocity are presented in Fig. 5.2 These temperature profiles are reported only for those cases that did not reach boiling. It can be observed that an increase in the inter-assembly gap sodium velocity leads to an overall decrease in the central assembly outlet bulk temperature, as expected. This dependence is considerably more prominent in the windowed (Fig. 5.2c) rather than the windowless (Fig. 5.2a) model, due to the mixing with the cooler sodium in the interassembly gap provided by the windows. This is further strengthened by comparing the differences in outlet bulk temperatures between the windowless and adiabatic (Fig. 5.2b) and windowed and adiabatic (Fig. 5.2 d) models. While these differences are relatively small at nominal operation (i.e. fC = 100%), they grow progressively larger for lower sodium velocities in the central assembly. This is consistent with

Table 4.3 Mesh convergence results for the outlet bulk temperature of the central assembly TC and gap TG for the high power case for the windowed web mesh. The results are presented for radial and axial refinements starting from the results obtained by a web mesh with: nR = 4, nS = 1, nG = 4 and a variable number of axial nodes depending on the axial region. Value [K]

Refinement

p[]

GCI1 [%]

GCI2 [%]

ACR [%]

Ex. val. [K]

Err. bnd. [K]

TC

864.56

Radial Axial

4.30 1.19

0.053 0.195

0.003 0.086

99.997 99.936

863.94 862.31

0.61 2.25

TG

712.64

Radial Axial

0.81 1.39

0.146 0.250

0.083 0.095

99.952 100.079

711.26 715.02

1.38 2.38

10

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

nominal operation the pressure difference between the gap and assemblies is negative, meaning that sodium outflows from the assemblies into the gap. With increasing gap velocity, the pressure drop through the axial length of the gap increases. The magnitude of the pressure difference between the gap and assemblies at the height of the windows thus decreases, limiting the outflow. If the pressure difference between gap and assembly was to be further increased and become positive by e.g. a decrease in the central assembly velocity, sodium inflow in the assembly is observed instead, leading to a reduction of the maximum temperatures with respect to the adiabatic model. With regard to the second observation, the non-monotonous behaviour of the maximum 40% is a consequence of temperature of the windowed model for fC the 3-Dimensionality of the sodium jet flowing into the central assembly from the windows. For fixed sodium velocity magnitudes at the inlet of the central assembly and the neighbouring assemblies, increasing the gap sodium velocity results in two effects: 1) the pressure difference between gap and central assembly increases, favouring sodium inflow into the central assembly; 2) the pressure in the vicinity of the windows in the central assembly decreases due to the increase of the velocity of the inflowing sodium jet. With regard to the second observation, this can be clearly observed in Fig. 5.4, as for a high sodium gap velocity (Fig. 5.3a), the velocity profile in the bulk of the central assembly deviates towards the windows, while this effect is considerably less prominent at low gap sodium velocity (Fig. 5.3b). From the perspective of variations in the maximum temperature in the assembly, these two effects are opposing. Sodium inflow in the assembly will favour the decrease in the maximum temperature. Conversely, a flow widening at the height of the window region will results in a lower axial centerline sodium velocity. Given that the maximum temperature is found to be on the assembly centerline at the height of the windows, this results in a increase of the maximum sodium temperature. This effect becomes more prominent at lower sodium velocities because of the increased velocity difference (and thus pressure difference) between the inflowing sodium jet from the windows and the central assembly sodium flow. With regard to the time-to-boil, it was found that the gap flow had no effect in appreciably delaying its onset in no case that reached boiling temperatures. This is due to the fact that at full power the heat 30% in the up process subsequent to a severe blockage (i.e. fC 20% in the windowed model) adiabatic and windowless models, fC is so rapid to be essentially adiabatic and in the order of magnitude of a few seconds. Thus, no results of interest are reported in this regard.

Fig. 4.3. Superficial velocity magnitude distribution at the height of the wrapper windows (black lines represent wrappers) for the high power case: a) for the coarsest mesh; b) for the finest mesh (i.e. number of nodes in all directions increased by a factor 4). The flow direction is dominantly oriented bottom-to-top.

60%, sodium inflow in the fact that, for sodium velocity fraction fC the central assembly from the windows is observed. This is observed even for a null inter-assembly gap sodium velocity, meaning that cooler sodium from the neighbouring assemblies leaks into central assembly undergoing blockage via the windows and inter-assembly gap. In particular, it was found that a central assembly sodium velocity fraction fC = 30% (i.e. 70% flow blockage) leads to boiling in the windowless model, but not in the windowed model due to the additional cooling windows provide. However, the dependence of the central assembly outlet bulk temperature on flow parameters does not necessarily represent the behaviour of the maximum temperature in the central assembly. As a matter of fact, it was found that the maximum central assembly temperature in the windowless model is essentially unaffected by gap flow (i.e. differences < 0.1 K ). This is not the case for the windowed model, as reported in Fig. 5.3. While the gap sodium velocity dependence of the maximum central assembly temperature is weaker than for the outlet bulk temperature (Fig. 5.3a), two observations can be made when comparing it against the adiabatic model (Fig. 5.3b): 1) at nominal operation ( fC = 100%) the maximum temperature in the windowed model is marginally higher than in the adiabatic model; 2) the behaviour of the maximum temperature with respect to gap velocity is not consistently monotonous. With regard to the first observation, this is due to the fact that at

5.2. Low power, low flow case A comparison of the central assembly outlet bulk temperatures, for varying degrees of central assembly sodium velocity and gap velocity are presented in Fig. 5.5. As before, these temperature profiles are reported only for those cases that did not reach boiling. The overall trends are the same as those observed for the full power case in Fig. 5.2. There is an overall increased in magnitude of the additional cooling power that the windowed (Fig. 5.5c) and windowless model (Fig. 5.5d) features provide with respect to the adiabatic case. The slope changes observed in the bulk temperature profiles for the windowless case (Fig. 5.5a and b) are associated with the laminarturbulent transition of the fluid flow in the inter-assembly gap. Outlet 30% because boiling is bulk temperatures are not reported for fC reached for both the windowless and windowed models. For the adia40%.With regard to the maxbatic model, boiling is reached for fC imum temperatures in the system and their dependence on flow parameters, results are presented in Fig. 5.6. At these low power and flow conditions, the inter-assembly gap is capable of sorting an effect on the maximum temperatures, unlike what was observed for the full power and flow case. This is due to the diffusive effects becoming more prominent as the importance of advective effects on heat transfer decreases at lower sodium velocities. 11

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

Fig. 5.1. Temperature distribution in the windowed model for fP = 5%, fN = 5%, fC = 40%, fG = 20% . On the right, the velocity streamlines are reported to highlight sodium leakage from the surrounding assemblies to the gap, and from the gap to the central assembly undergoing partial blockage. The wrappers are highlighted in white.

Fig. 5.2. Central assembly outlet temperatures for the: a) windowless model; c) windowed model. Difference between the central assembly temperatures between: b) the windowless and adiabatic models; d) the windowed and adiabatic model. The temperatures are reported against the gap sodium velocity fraction fG for different central assembly sodium velocities fC . 12

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

Fig. 5.3. High power, high flow case: a) maximum temperature in the central assembly of the windowed model; b) difference between the maximum temperatures obtained by the windowed and adiabatic models. These are presented against inter-assembly gap sodium velocity fraction fG at various central assembly sodium velocity fractions fC .

6. Conclusions and future developments The present work dealt with the further expansion of the GeN-Foam porous-medium based coarse-mesh thermal-hydraulic solver in order to characterise some features considered in the latest ESFR-SMART core design, namely novel assembly wrapper windows and the inter-assembly gap, and to assess their relative performance. A novel wrapper heat transfer model was developed and verified against analytical solutions. The developed heat transfer model neglects the thermal inertia of the wrapper, which is acceptable for thin baffles made of a relatively low specific heat capacity material so that the absolute heat capacity ratio between the sodium and wrapper is small (such was the case in our investigation). The localized pressure losses at the windows were modelled as those that would be expected from a 90-degree miter bend of the same hydraulic diameter. In order to characterise the thermalhydraulic impact of the investigated features, three different geometric models were produced: 1) an axial section of a single assembly spanning only the axial fuel region with thermally insulating wrappers (referred to as “adiabatic” model); 2) an axial section of a bundle of seven assemblies with thermally conductive wrappers spanning only the axial fuel region (referred to as “windowless” model); 3) an axial section of a bundle of seven assemblies with thermally conductive wrappers and windows spanning all the axial regions found between the fuel region and the assembly outlet (referred to as “windowed” model). The results consisted in the coolant bulk temperatures at the assembly and gap outlets and the maximum coolant temperature in the system. These were obtained for varying degrees of assembly power, assembly sodium velocity (to model different core flow levels and/or flow blockage in the central assembly) and inter-assembly gap sodium velocity. A mesh convergence study was performed to assess the most suitable coarse mesh for each of the three models and to establish the magnitude of the expected uncertainty bands. A coarse mesh has a limited room for improvement in terms of mesh independence of the results when compared to traditional fine mesh approaches. However, as long as it can be shown that the results are at least in the asymptotic convergence range, results extrapolation at null mesh spacing is possible, as well as uncertainty band estimation. For two representative sets of flow and power conditions, the convergence properties of the outlet bulk temperatures were assessed for the windowed and windowless model on a variety of tentative meshes. The chosen meshes satisfied the requirements of providing outlet bulk temperatures in the asymptotic convergence range. The estimated uncertainty bands on the outlet bulk temperatures varied in the 1 K 8 K range depending on the case. The parametric investigation at two different power and flow levels was performed on each of the three models. In particular, these

Fig. 5.4. Comparison between the full power, full flow velocity streamlines at fC = 30% in the windowed model for two different inter-assembly gap velocities: a) fG = 100% , b) fG = 20% . The view consists of a 2-D slice parallel to the main flow direction that encompasses the central assembly, the inter-assembly gap and part of the neighbouring assemblies. The position of the assembly wrappers and thus the windows is highlighted in black.

With regard to the time-to-boil, both the windowless and windowed models were found to perform comparably in term of delaying the onset of boiling with respect to the adiabatic case. For this reason, only timeto-boil results relative to the windowed model and its comparison against the adiabatic model are presented in Fig. 5.7. The overall boiling delay effects grow smaller in magnitude for lower central assembly sodium velocities (i.e. higher blockages) as expected. On the other hand, inter-assembly gap flow can provide an additional margin to boiling delay. Nonetheless, the magnitude of the delay is of negligible engineering importance.

13

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

Fig. 5.5. Central assembly outlet temperatures in the low power, low flow case for the: a) windowless model; c) windowed model. Difference between the central assembly temperatures between: b) the windowless and adiabatic models; d) the windowed and adiabatic model. The temperatures are reported against the gap sodium velocity fraction fG for different central assembly sodium velocities fC .

states consisted of: 1) nominal power and flow conditions; 2) low power and flow conditions compatible with core state after a protected transient and onset of natural circulation. For each state, a parametric investigation was performed with respect to inter-assembly gap sodium velocity and central assembly sodium velocity, to model flow blockage. The most important results are hereby summarised. 1) At both nominal and low (i.e 5% of nominal) flow conditions, the axial structure of the assembly shapes the pressure drop profiles in such a way that sodium outflow from the assembly is observed. This results in a small assembly outlet bulk temperature increase 5 K depending on the inter-assembly gap sodium velocity. 2) For any investigated central assembly sodium velocity below nominal values, sodium inflow in the assembly through the wrapper windows is observed. This was always found to reduce both the assembly outlet temperatures and the maximum assembly temperature by appreciable quantities with respect to the thermally insulated, windowless wrappers. 3) The change in maximum assembly coolant temperature with respect to inter-assembly gap sodium velocity is not always monotonous in the presence of wrapper windows. This was found to be due specifically to 3-D effects of the inflowing sodium jet from the wrapper windows. 4) Inter-assembly gap flow in the absences of wrapper windows was found to be less effective than the windowed wrapper in lowering outlet bulk temperatures for varying degrees of inter-assembly gap sodium velocity. Furthermore, it was found to have no effect at all on the maximum assembly temperatures at nominal flow conditions due to the dominance of advective effects compared to diffusive effects. 5) It was found that inter-assembly gap flow, with or without windows, is capable of preventing boiling for

some flow conditions at nominal power levels when compared to the thermally insulated wrapper case. It is recalled that the main purpose of the assembly windows is to limit the void presence to the upper fuel region in case of boiling. A full assessment of the performance of the wrapper windows thus requires two-phase flow modelling and coupling with the neutronics sub-solver of GeN-Foam. The expansion of the currently single-phase coarse-mesh solver to two-phase flow analysis is currently being undertaken and will enable a more detailed investigation of the wrapper windows. On the single-phase side, future work consists in improving the physical models used for the description of momentum transfer, heat transfer and turbulence. Such improvements will rely both on state-of-the-art correlations for pressure drops (e.g. (Chen et al., 2014)) and fluid-wrapper wall heat transfer coefficients. Improvements on the turbulence modelling side will mainly rely on the inclusion of experimental data to fine-tune the turbulent thermal diffusivity to account for the presence of wrapping wire in a more accurate manner. Efforts concerning the verification of the solver capabilities against a modified version of the TRACE system code (Chenu et al., 2009) have been initiated and will be presented together with future planned validation efforts. With regard to solver validation efforts, the Laboratory for Reactor Physics and Systems Behaviour is currently taking part in an IAEA coordinated research project related to the benchmark analysis of past experiments conducted at the Fast Flux Test Facility (FFTF) experimental sodium fast reactor. In particular, the investigation is centered on the Loss Of Flow Without SCRAM (LOFWOS) (IAEA, 2019) experiments that represent an excellent opportunity for multi-physics code validation and code-to-code comparison efforts. 14

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

Fig. 5.6. Maximum temperature in the low power, low flow case for the: a) windowless model; c) windowed model. Differences of the maximum temperature between: b) the windowless and the adiabatic models; d) the windowed and the adiabatic models. These are presented against inter-assembly gap sodium velocity fraction fG at various central assembly sodium velocity fractions fC .

Fig. 5.7. Low power, low flow case: a) time-to-boil of the windowed model; b) difference in time-to-boil between the windowed and adiabatic models. These are presented against the inter-assembly gap sodium velocity fraction fG for various central assembly sodium velocity fractions fC .

15

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

Declaration of Competing Interest

Acknowledgements

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.

The authors would like to acknowledge the cooperation with the Horizon 2020 ESFR-SMART project (Mikityuk et al., 2017). Geometries and meshes have been generated with ad hoc developed scripts for the Trelis meshing software by csimsoft™ (2019).

Appendix A. EFSR-SMART data The latest design specifications of the ESFR-SMART (Mikityuk et al., 2017) core that have been employed for the creation of the computational domains can be found in Table A.1. The material properties employed in the simulations (Fink and Leibowitz, 1996) within the scope of this work are reported in Table A.2. Table A.1 Reference data of the latest ESFR-SMART core Quantity

Value

Reactor thermal power [MWth] Number of inner fuel assemblies Number of outer fuel assemblies Number of pins per fuel assembly Core outlet temperature [° C] Core inlet temperature [° C]

3600 216 288 271 395 545 19000

Mass flow rate through active core

kg s

Fuel material Fuel enrichment (%wt) Cladding and wrapper material Active fuel region height (inner core) [mm] Active fuel region height (outer core) [mm] Assembly lattice pitch [mm] Wrapper outer flat-to-flat [mm] Wrapper thickness [mm] Inter-assembly gap thickness [mm] Fuel pellet inner diameter [mm] Fuel pellet outer diameter [mm] Cladding inner diameter [mm] Cladding outer diameter [mm]

(U, Pu) O2 17.99 EM10 Steel 750.00 950.00 209.85 205.35 4.50 4.50 1.56 4.68 4.84 5.34

Table A.2 Reference material properties. Material

Quantity

Sodium

841.06

kg m3

J kg K

cp (T )

W mK

k (T )

µ (T ) [Pa s] Fuel

k Cladding

J kg K W mK

k

J kg K W mK

h gap

4.472·10

3

2.457·10 5T + 6.427·10 8T2

·10 17T 5 7.603·10 21T 6 124.670 1.138·10 1T + 5.523·10 5T2

250 2 7500

kg m3

cp

6.322·10 1T + 2.254·10 4T2 + 7.515·10 8T3

1583.000

10542

kg m3

cp

Value

W m2 K

500 20 3000

16

9.486·10

11T3

1.184·10 8T3

+ 8.150·10

14T 4

3.820

Nuclear Engineering and Design 355 (2019) 110291

S. Radman, et al.

References

https://doi.org/10.1016/j.anucene.2010.06.010. Boure, J., Bergles, A., Tong, L., 1973. Review of two-phase flow instability. Nucl. Eng. Des. 25 (2), 165–192. https://doi.org/10.1016/0029-5493(73)90043-5. Kottowski, H., Savatteri, C., 1984. Fundamentals of liquid metal boiling thermohydraulics. Nucl. Eng. Des. 82 (2), 281–304. https://doi.org/10.1016/00295493(84)90216-4. Yue, N., Zhang, D., Chen, J., Song, P., Wang, X., Wang, S., Qiu, S., Su, G., Zhang, Y., 2018. The development and validation of the inter-wrapper flow model in sodium-cooled fast reactors. Prog. Nucl. Energy 108, 54–65. https://doi.org/10.1016/j.pnucene. 2018.05.007. Droin, J.-B., Marie, N., Bachrata, A., Bertrand, F., Merle, E., Seiler, J.-M., 2017. Physical tool for unprotected loss of flow transient simulations in a sodium fast reactor. Ann. Nucl. Energy 106, 195–210. https://doi.org/10.1016/j.anucene.2017.03.035. Hu, R., Yu, Y., 2016. A computationally efficient method for full-core conjugate heat transfer modeling of sodium fast reactors. Nucl. Eng. Des. 308, 182–193. https://doi. org/10.1016/j.nucengdes.2016.08.018. Radman, S., Fiorina, C., Pautz, A., 2017. Investigation of partial coolant flow blockage in a sodium fast reactor assembly with coarse-mesh methodologies. In: Proceedings of 26th International Conference New Energy for New Europe. Bertini, H.W., 1980. Descriptions of selected accidents that have occurred at nuclear reactor facilities. doi:https://doi.org/10.2172/5591584. Vafai, K., 2005. Handbook of Porous Media, second ed. CRC Press. Whitaker, S., 1985. A simple geometrical derivation of the spatial averaging theorem. Chem. Eng. Educ. 19 (1). Chandesris, M., Serre, G., Sagaut, P., 2006. A macroscopic turbulence model for flow in porous media suited for channel, pipe and rod bundle flows. Int. J. Heat Mass Transfer 49 (15), 2739–2750. https://doi.org/10.1016/j.ijheatmasstransfer.2005.12. 013. Nakayama, A., Kuwahara, F., 1999. A macroscopic turbulence model for flow in a porous medium. J. Fluids Eng. 121 (2), 427–433. https://doi.org/10.1115/1.2822227. Launder, B.E., Sharma, B., 1974. Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc. Lett. Heat Mass Transfer 1 (2), 131–137. ANSYS FLUENT 12.0 user manual. URL:http://www.afs.enea.it/project/neptunius/docs/ fluent/html/ug/main_pre.htm. Khan, E., Rohsenow, W., Sonin, A., Todreas, N., 1975. A porous body model for predicting temperature distribution in wire-wrapped fuel rod assemblies. Nucl. Eng. Des. 35 (1), 1–12. https://doi.org/10.1016/0029-5493(75)90076-X.. URL:http://www. sciencedirect.com/science/article/pii/002954937590076X. Holzmann, T., 2018. Mathematics, numerics, derivations and OpenFOAM. URL:https:// holzmann-cfd.de/publications/mathematics-numerics-derivations-and-openfoam/ item/17-mathematics-numerics-derivations-and-openfoam. Rehme, K., 1973. Pressure drop correlations for fuel element spacers. Nucl. Technol. 17, 12–23. Gunter, A., Shaw, W., 1945. A general correlation of friction factors of various types of surfaces in cross flow. ASME Trans. 67, 643–660. Mikityuk, K., 2009. Heat transfer to liquid metal: review of data and correlations for tube bundles. Nucl. Eng. Des. 239 (4), 680–687. https://doi.org/10.1016/j.nucengdes. 2008.12.014. Yu, Y., Merzari, E., Obabko, A., Thomas, J., 2015. A porous medium model for predicting the duct wall temperature of sodium fast reactor fuel assembly. Nucl. Eng. Des. 295, 48–58. https://doi.org/10.1016/j.nucengdes.2015.09.020.. URL:http://www. sciencedirect.com/science/article/pii/S0029549315004276. Al-Tameemi, W., Ricco, P., 1975. Pressure-loss coefficient of 90 deg sharp-angled miter elbows. J. Fluids Eng. 140, 1–12. Roache, P.J., 1994. A method for uniform reporting of grid refinement studies. J. Fluids Eng. 116, 405–413. https://doi.org/10.1115/1.2910291. Sidi, A., 1997. A complete convergence and stability theory for a generalized richardson extrapolation process. SIAM J. Numer. Anal. 34 (5), 1761–1778. https://doi.org/10. 1137/S0036142994278589. Chen, S., Todreas, N., Nguyen, N., 2014. Evaluation of existing correlations for the prediction of pressure drop in wire-wrapped hexagonal array pin bundles. Nucl. Eng. Des. 267, 109–131. https://doi.org/10.1016/j.nucengdes.2013.12.003.. URL:http:// www.sciencedirect.com/science/article/pii/S002954931300678X. Chenu, A., Mikityuk, K., Chawla, R., 2009. TRACE simulation of sodium boiling in pin bundle experiments under loss-of-flow conditions. Nucl. Eng. Des. 239 (11), 2417–2429. https://doi.org/10.1016/j.nucengdes.2009.07.015. IAEA. Coordinated research project on FFTF LOFWOS tests. URL:https://www.iaea.org/ newscenter/news/new-crp-benchmark-analysis-of-fftf-loss-of-flow-without-scramtest-i32011. csimsoft™. Trelis. URL:https://www.csimsoft.com/trelis.jsp. Fink, J., Leibowitz, L., 1996. A consistent assessment of the thermophysical properties of sodium. High Temp. Mater. Sci. 35 (1). URL:https://www.osti.gov/biblio/251097.

OECD-NEA, 2014. Technology roadmap update for Generation IV nuclear energy systems. URL:https://www.gen-4.org/gif/upload/docs/application/pdf/2014-03/gif-tru2014. pdf. Tenchine, D., 2010. Some thermal hydraulic challenges in sodium cooled fast reactors. Nucl. Eng. Des. 240 (5), 1195–1217. https://doi.org/10.1016/j.nucengdes.2010.01. 006. CEA, 1978. La centrale à neutrons rapides de Creys-Malville (Super-Phénix). Bull. Inform. Scientifiques Tech. 227, 139. Lefvre, J., Mitchell, C., Hubert, G., 1996. European fast reactor design. Nucl. Eng. Des. 162 (2), 133–143. https://doi.org/10.1016/0029-5493(95)01141-2. Fiorini, G., Vasile, A., 2011. European commission – 7th framework programme: the collaborative project on european sodium fast reactor (CP ESFR). Nucl. Eng. Des. 241 (9), 3461–3469. https://doi.org/10.1016/j.nucengdes.2011.01.05. Seventh European Commission conference on Euratom research and training in reactor systems (Fission Safety 2009). Mikityuk, K., Girardi, E., Krepel, J., Bubelis, E., Fridman, E., Rineiski, A., Girault, N., Payot, F., Buligins, L., Gerbeth, G., Chauvin, N., Latge, C., Garnier, J.-C., 2017. ESFRSMART: new Horizon-2020 project on SFR safety. In: Proceedings of International Conference on Fast Reactors and Related Fuel Cycles: Next Generation Nuclear Systems for Sustainable Development, no. IAEA-CN245-450. OpenCFD Ltd., OpenFOAM. URL:https://www.openfoam.com/about/. Weller, H.G., Tabor, G., 1998. A tensorial approach to computational 10 continuum mechanics using object-oriented techniques. Comput. Phys. 12, 620. https://doi.org/ 10.1063/1.168744. Jasak, H., 2009. OpenFOAM: open source CFD in research and industry. Int. J. Naval Arch. Ocean Eng. 1 (2), 89–94. https://doi.org/10.2478/IJNAOE-2013-0011. Fiorina, C., Clifford, I., Aufiero, M., Mikityuk, K., 2015. GeN-Foam: a novel OpenFOAM based multi-physics solver for 2D/3D transient analysis of nuclear reactors. Nucl. Eng. Des. 294, 24–37. https://doi.org/10.1016/j.nucengdes.2015.05.035. Wahnon, S.S.P., Ammirabile, L., Kloosterman, J., Lathouwers, D., 2018. Multi-physics models for design basis accident analysis of sodium fast reactors. Part i: validation of three-dimensional TRACE thermal-hydraulics model using Phénix end-of-life experiments. Nucl. Eng. Des. 331, 331–341. https://doi.org/10.1016/j.nucengdes. 2018.02.038. Bottoni, M., Dorr, B., Homann, C., Struwe, D., 1987. State of development of the computer programme BACCHUS-3D/TP for the description of transient two-phase flow conditions in LMFBR fuel pin bundles. Nucl. Eng. Des. 100 (3), 321–349. https://doi.org/ 10.1016/0029-5493(87)90084-7. Chavez, V.J., Imke, U., Sanchez-Espinoza, V., 2018. TWOPORFLOW: a two-phase flow porous media code, main features and validation with BWR-relevant bundle experiments. Nucl. Eng. Des. 338, 181–188. https://doi.org/10.1016/j.nucengdes.2018.08. 009. Jeong, J., Yoon, H., Park, I., Cho, H., 2010. The CUPID code development and assessment strategy. Nucl. Eng. Technol. 42 (6), 636–655. https://doi.org/10.5516/NET.2010. 42.6.636. Fiorina, C., Kerkar, N., Mikityuk, K., Rubiolo, P., Pautz, A., 2016. Development and verification of the neutron diffusion solver for the GeN-Foam multi-physics platform. Ann. Nucl. Energy 96, 212–222. https://doi.org/10.1016/j.anucene.2016.05.023. Fiorina, C., Hursin, M., Pautz, A., 2017. Extension of the GeN-Foam neutronic solver to SP3 analysis and application to the CROCUS experimental reactor. Ann. Nucl. Energy 101, 419–428. https://doi.org/10.1016/j.anucene.2016.11.042. de Oliveira, R.G., Mikityuk, K., 2018. Analytical solutions to a coupled fluid dynamics and neutron transport problem with application to GeN-Foam verification. Ann. Nucl. Energy 121, 446–451. https://doi.org/10.1016/j.anucene.2018.07.036. Radman, S., Fiorina, C., Mikityuk, K., Pautz, A., 2018. A simplified numerical benchmark for pool-type sodium fast reactors. In: 26th International Conference on Nuclear Engineering. American Society of Mechanical Engineers. https://doi.org/10.1115/ ICONE26-82260. pp. V004T15A018–V004T15A018. Bortot, S., Mikityuk, K., Panadero, A., Pelloni, S., Alvarez-Velarde, F., Lopez, D., Fridman, E., Cruzado, I., Herranz, N., Ponomarev, A., 2015. European benchmark on the ASTRID-like low-void-effect core characterization: neutronic parameters and safety coefficients. In: International Congress on Advances in Nuclear Power Plants, . URL:http://inis.iaea.org/search/search.aspx?orig_q=RN:49004848. Mikityuk, K., Chenu, A., Sun, K., 2011. European patent: a wrapper tube for a fuel subassembly of a nuclear reactor core and method for protecting fuel against overheating in case of coolant boiling. URL:https://patentscope.wipo.int/search/en/detail.jsf? docId=WO2013098079. Sun, K., Krepel, J., Mikityuk, K., Pelloni, S., Chawla, R., 2011. Void reactivity decomposition for the sodium-cooled fast reactor in equilibrium fuel cycle. Ann. Nucl. Energy 38 (7), 1645–1657. https://doi.org/10.1016/j.anucene.2011.02.018. Tommasi, J., Archier, P., Ruggieri, J., 2010. Validation of the sodium void reactivity effect prediction using JEFF-3.1 nuclear data. Ann. Nucl. Energy 37 (11), 1534–1553.

17