# STH-CFD Codes Coupled Calculations Applied to HLM Loop and Pool Systems

Authors: M. Angelucci; D. Martelli; G. Barone; I. Di Piazza; N. Forgione
Journal: Science and Technology of Nuclear Installations
 (2017)

---

## Abstract

This work describes the coupling methodology between a modified version of RELAP5/Mod3.3 and ANSYS Fluent CFD code developed at the University of Pisa. The described coupling procedure can be classified as “two-way,” nonoverlapping, “online” coupling. In this work, a semi-implicit numerical scheme has been implemented, giving greater stability to the simulations. A MATLAB script manages both the codes, oversees the reading and writing of the boundary conditions at the interfaces, and handles the exchange of data. A new tool was used to control the Fluent session, allowing a reduction of the time required for the exchange of data. The coupling tool was used to simulate a loop system (NACIE facility) and a pool system (CIRCE facility), both working with Lead Bismuth Eutectic and located at ENEA Brasimone Research Centre. Some modifications in the coupling procedure turned out to be necessary to apply the methodology in the pool system. In this paper, the comparison between the obtained coupled numerical results and the experimental data is presented. The good agreement between experiments and calculations evinces the capability of the coupled calculation to model correctly the involved phenomena.


---

## Body

## 1. Introduction

Thermal-hydraulic analyses are a fundamental issue in the development, design, and licensing of nuclear power plants (NPPs). The investigation of the plant performance during accidental conditions has always been one of the main concerns of nuclear safety. System Thermal-Hydraulic (STH) codes, such as RELAP5, CATHARE, ATHLET, and TRACE, are widely used in the assessment and improvement of the safety aspects of existing and new NPPs. The Best-Estimate (BE) analysis with STH codes allows the simulation of accidental transients as realistically as possible and is used to assess the values of safety margins. BE system codes are currently the only ones accepted for licensing purposes, provided that they are properly qualified and used together with a methodology to evaluate uncertainties.STH codes are based on two-phase flow, one-dimensional (1D) equations, where mass, momentum, and energy conservation equations are solved for both the liquid and vapor phases separately. Several constitutive equations and empirical correlations are required to close the set of balance equations [1]. The optimization and the improvement of the tools used for thermal-hydraulic safety analysis are an ongoing process and a permanent objective of the research efforts. One of the major limits of the STH codes is related to its one-dimensional approach, which prevents the study of more complex 3D phenomena. In this context, the use of computational fluid dynamic (CFD) codes has increased in the last fifteen years. CFD codes provide more detailed information, since the discretization of the domain is highly refined. They also include models for simulating turbulence, heat transfer, multiphase flows, and chemical reactions. Nevertheless, the complexity of the used models and the high spatial resolution results in higher computational efforts compared with those required by the STH codes when applied to an equivalent domain.In this scenario, coupling techniques between STH and CFD codes for thermal-hydraulic analysis purposes are emerging worldwide. The main advantage of this methodology is the possibility of modelling multiscale phenomena, at different level of detail. The CFD code should be used to analyse a smaller part of the domain where 3D effects are significant and/or detailed flow information is demanded. The STH code should be applied to the portions of the system characterized by 1D components (i.e., pipes) and to model multiphase and/or multicomponent flow.First developments and applications of coupling methodology can be found in Korea, in the works of Jeong et al. [2, 3], and in USA, in the works of Aumiller et al. [4, 5]. More recent works on coupled calculation applied to water systems can be found in literature [6–9]. Moreover, coupling methodologies found wide application in GEN-IV systems thermal-hydraulic analysis, as in [10] where coupled simulation was performed on a gas-cooled reactor. In Europe, the French Atomic Energy Commission (CEA) developed a coupling tool between the 3D CFD code TRIO_U with the BE-STH code CATHARE, with applications to PHENIX Natural Circulation Test [11]. GRS Society for Plant and Reactor Safety (Germany) worked on the coupling of the 1D System Code ATHLET with the 3D CFD Software Package ANSYS CFX [12], while KTH Royal Institute of Technology developed a coupling methodology between RELAP5 and the CFD code STAR-CCM+. In the context of the 7th Framework EU Program THINS project, GRS and KTH applied the developed tools to the Heavy Liquid Metal (HLM) system TALL-3D. The comparison of the two coupling approaches is reported in [13]. At the Belgian Nuclear Research Centre (SCK·CEN), a numerical algorithm of a STH/CFD coupling method was developed for multiscale transient simulations of pool-type reactors [14]. The application of STH codes to GEN-IV system required the implementation of thermophysical properties and correlations for new kind of fluids (e.g., liquid metals). Balestra et al. at the University of Rome, in collaboration with the Idaho National Laboratory (INL), developed a tool for the generation of new properties binary file needed by the code RELAP5-3D [15] and implemented thermophysical properties for LBE and Lead recommended in [16].At the University of Pisa (UniPi), the relationship used by the RELAP5/Mod3.3 to generate the table of lead, LBE, lead-lithium, and sodium properties were modified [17] in agreement with the ones suggested in [18] and then taken as the reference in [16]. Further, specific convective heat transfer correlations for LMs have been implemented inside the code. This modified version of RELAP5/Mod3.3 was adopted at the University of Pisa in order to develop an in-house coupling methodology between the STH code and Fluent commercial CFD code since 2012. The first release of the coupling tool used an explicit scheme for the advancement in time. This tool was applied to the NACIE facility, with good agreement of the obtained results against experimental data [19, 20]. In the present work, the improvements of the coupling methodology developed at UniPi, such as the semi-implicit numerical scheme, are presented in Section 2. The coupling tool has been applied to NACIE and CIRCE, respectively, a loop and a pool facility designed and built at the Italian National Agency for New Technologies, Energy and Sustainable Economic Development (ENEA) in Brasimone RC, both using a HLM as working fluid. The computational domain of the two facilities and the main results of the coupled calculation are presented, respectively, in Sections 3 and 4. The main conclusions of this work are summarized in Section 5.


## 2. Coupling Procedure

The coupling approach consists in a “nonoverlapping,” “two-way” coupling scheme. In a “nonoverlapping” strategy the overall domain is divided into some regions modelled using the CFD approach(A) and other regions modelled using a system code (B), similar to the simplified scheme reported in Figure 1. From the geometrical domain subdivision, one or more “coupling” interfaces (Π) are conveniently selected; here the thermofluid-dynamics data are transferred from the CFD code to the STH-code-region and vice versa (“two-way coupling”). Several and different physical parameters (e.g., mass flow, velocity, pressure, and temperature) can be passed at each interface Π. The number and the type of exchanged quantities are not necessary the same for each interface (referring to Figure 1  p not necessarily equal to q-p and/or r-q).Figure 1
Simplified sketch of a “nonoverlapping,” “two-ways” coupled scheme.If the set of CFD and STH balance equations were merged into a single system and solved by a unique code, a so-called “monolithic solution” would be implemented. Yet, this method requires major modifications to the source codes (if available). Alternatively, a “partitioned solution” can be adopted, where each code solves independently its own set of equations and a coupling interface is required to exchange thermal-hydraulic variables at the interfaces. In this work, the “partitioned solution” approach was chosen and a third software, MATLAB®, was selected to manage the coupling interfaces. An appropriate MATLAB script handles the codes synchronization, drives the execution of both the solvers, enables the exchange, and processes the data when required.In the previous work of Martelli et al. [20], an explicit numerical scheme was implemented. In Figure 2, the logic of the explicit (a) scheme compared to the semi-implicit (b) scheme is illustrated. Regarding the explicit scheme, at the beginning of the run (n=1), the Fluent (master code) advances first from the initialization timet0 = 0 up to the end of the first time step (Δt). Suitable CFD results are used as BCs for a RELAP5 simulation. Accordingly, the RELAP5 (slave code) computation follows and its results are temporarily saved in MATLAB, in order to be used as Fluent BC in the following time step. At the end of both the runs (Fluent and RELAP5), the coupled calculation transient time advances (n=n+1). The main drawback of the explicit method is that the conservation of mass, momentum, and energy at the interfaces is not necessarily guaranteed and very small time step could be mandatory to avoid numerical instabilities. The time step is usually chosen considering the involved physical phenomena, always respecting the Courant–Friedrichs–Lewy (CFL) condition for the stability. The chosen value is used as time step in the CFD calculation and as duration of each RELAP5 transient calculation. To reduce the computational cost and improve the numerical stability, a semi-implicit numerical scheme was developed and implemented in the present work. Since the two codes solve independently their own set of equations (“partitioned solution”), it is not possible to implement a fully implicit numerical scheme, which needs the solution of a single system of discretized equations. Nevertheless, it is common practice to refer to this semi-implicit method as “implicit,” as it could be done later on. In the semi-implicit scheme (Figure 2(b)), a generic external iteration (n) begins with the time step advancement Δt. The choice of the time step is mainly guided by the CFD portion of the calculation, considering the physical phenomena involved. The same time step is used in RELAP5 as total duration of each transient calculation. An inner iteration (j=1) also starts with the calculation of the master code Fluent. The CFD results are set as BCs for the STH code RELAP5, which runs for the same time step, and its results are read by MATLAB and stored in a “txt” file. At this point, the convergence is checked by comparing the updated value of all thermal-hydraulic quantities xi at each interface Π with the results of the previous iteration stored in MATLAB. If the convergence criterion is satisfied, a new external iteration n+1 begins, and the updated RELAP5 results are transmitted to a new Fluent calculation (from n to n+1). Otherwise, a new inner iteration is performed (j=j+1) and the last RELAP5 results are transmitted again to the Fluent calculation related to the time step n. The convergence parameter εi is usually the same for each type of physical parameter (mass flow, velocity, temperature, pressure, etc.) and is kept constant during the whole simulation. This approach assures the continuity at the physical interfaces and the conservation of the main physical parameters (mass, momentum, and energy). In fact, each code is verified independently, whereas the succession of inner iterations ensures the coherence between the information transferred from the CFD outlet boundary to the RELAP5 inlet boundary and the information transferred from the RELAP5 outlet boundary to the RELAP5 inlet boundary.Figure 2
Logic of the explicit (a) and semi-implicit (b) numerical schemes.
(a)
(b)A new logic, developed to manage the Fluent code from MATLAB, is also illustrated in the present work. The basic idea of the new methodology is to execute the Fluent code directly from MATLAB at the beginning of the simulation in Server mode (“Fluent as a Server”) on the Local PC used. The use of CORBA add-ons for Fluent [21] allows the establishment of an interface between MATLAB and “Fluent as a Server” so that it is possible to manage the Fluent session through MATLAB Text User Interface (TUI) specific commands. This method brought about a great advantage in terms of global computational time with respect to the method used before (about 4 times faster), avoiding, furthermore, the use of a parallelized User Define Function (UDF) and journal files to exchange data between Fluent and RELAP5.Following the scheme of Figure3, after the CFD computation, the MATLAB script is in charge of reading the results necessary as BCs for RELAP5. The above-mentioned data are written in the input file for RELAP5 calculation, which is then launched. Similar to what is described before, the results from RELAP5 computation are read through the MATLAB script and the information needed to be exchanged with Fluent is extracted. Then, the MATLAB process checks the convergence and, consequently, a new inner or external iteration occurs. RELAP5 and Fluent files are saved and overwritten at the end of each inner iteration, whereas converged calculations and main results are saved at the end of each time step (outside the inner iterations).Figure 3
Coupled code management through the MATLAB software.


## 3. Application to Loop Facility

### 3.1. NACIE Facility

NACIE [22] is a Lead Bismuth Eutectic (LBE) loop-type facility. It consists of a rectangular loop made of two vertical stainless steel pipes (2.5′′, Sch. 40) 7.5 m long, acting as riser and downcomer, connected with two horizontal pipes 1 m long. The layout of the facility is shown in Figure 4. A heat source is installed in the bottom part of the riser and consists of two electrical heated rods with nominal thermal power of about 43 kW. A water/LBE “tube in tube” counterflow heat exchanger (HX) is placed on the upper part of the downcomer. The secondary side is fed by water at low pressure. The facility can work in both natural and gas enhanced circulation conditions. Indeed, argon can be injected downstream the heat source, through a 10 mm diameter pipe, housed inside the riser from the expansion vessel.Figure 4
Isometric view of NACIE primary loop.The computational domain for RELAP5/Fluent coupled calculations of the NACIE facility, already described and used in previous work [20], is depicted in Figure 5. It is divided into two nonoverlapping regions: the CFD code was used to simulate the Fuel Pin Simulator (FPS) representing the heat source, whereas the remaining parts of the loop and the HX (both primary and secondary sides) were modelled with the STH code RELAP5. Referring to Figure 5(a),Pipe-130 represents the riser, whileTDV-400, TDJ-405, andBranch-125, at bottom ofPipe-130, simulate the injection of the argon which promote the circulation of the LBE. The expansion vessel (fromPipe-146 toPipe-156) allows the separation between argon and LBE, so that the LBE can flow through the horizontal pipe and then through the HX(Pipe-180) and the downcomer (Pipe-200 andPipe-206). The secondary side water system was modelled by means ofTDV-500,TDJ-505, the HX secondary side annular zone(Annulus-510), andTDV-520.Figure 5
NACIE loop computational domain (a); codes exchanged data (b).
(a)
(b)The FPS was modelled by the Fluent code with both a simplified 2D axial-symmetric domain and a 3D model, shown, respectively, in Figures6(a) and 6(b). The 2D geometrical model was discretized by a structured mesh composed of 7668 rectangular cells uniformly distributed in both the axial and radial coordinates. To model the FPS form losses, mainly due to spacer grids, a constant pressure drop coefficient was assumed in the 2D domain. For this purpose, five distinct interior faces were set as “porous jump,” with each one characterized by an equivalent constant local pressure drop coefficient of 0.7. The 3D geometrical domain was modelled using a symmetry plane, as shown in Figure 6(b). The electrical pins were not included in the domain and heat flux boundary condition (BC) was used at the pin walls. The domain was discretized using 141045 hexahedral elements with refinement close to pins wall, near the inlet and outlet sections. One interior face was set as a porous jump and an equivalent constant coefficient of concentrated pressure drop equal to 0.5 was set in order to simulate the pressure drop due to the spacer grid not simulated in the CFD model.Figure 6
2D (a) and 3D (b) CFD domains of the Fuel Pin Simulator (FPS).
(a)
(b)Concerning the coupled calculation, the exchange of data takes place at the inlet and outlet of the FPS, following the scheme displayed in Figure5(b). Fluent gives the mass flow m˙2 and temperature T2, at the outlet of the FPS, to RELAP5 inTDJ-115 andTDV-112; the pressure p1, obtained from the inlet section of the CFD domain, is set inTDV-110 in the RELAP5 domain. Similarly, the CFD BCs at inlet section, that is, mass flow m˙1 and temperature T1, are obtained from the values computed from RELAP5 in the last volume ofPipe-100, while the outlet surface BC pressure p2 is obtained from the value computed in the first volume ofPipe-120.


### 3.2. Test Matrix

The application to NACIE facility of RELAP5/CFD coupled calculations concerned an isothermal test in gas enhanced circulation (Test 206). During the experimental test, the FPS was switched off and the HX was empty and nonactive. The LBE temperature remained in the range of 200–250°C. The experimental test matrix is shown in Table1, which reports the adopted boundary conditions. Test 206 started with LBE at rest in the facility (0 kg/s mass flow rate); then the gas mass flow controller was switched on to a nominal value of 2 Nl/min (Step 1). The gas flow rate was kept constant for the time needed to have a steady state statistically significant. After that, the other steps came in succession, increasing or decreasing the Ar flow rate with the order delineated in Table 1.Table 1
Boundary conditions of Test 206 in NACIE facility.
T
L
B
E
FPS power
HX power
Nominal Ar flow rate [Nl/min]
200–250°C
0 kW
0 kW
Step 1
2.0
Step 7
8.0
Step 2
4.0
Step 8
6.0
Step 3
5.0
Step 9
5.0
Step 4
6.0
Step 10
4.0
Step 5
8.0
Step 11
2.0
Step 6
10.0


### 3.3. Results

The coupled calculation started with the facility at rest and no gas injected in the riser. After the injection of argon was activated (flow rate to 2 Nl/min), the LBE mass flow rate increased up to a value of about 7.7 kg/s, as shown in the first step of Figure7(a). Increasing the argon flow rate (4-5-6-8-10 Nl/min), the LBE mass flow rate increased up to a maximum of about 14 kg/s. In the second part of the test, gas injection was decreased symmetrically with respect to the increasing ramp and the LBE flow rate followed the same trend. As it can be noticed in Figure 7(a), calculated LBE mass flow rate tended to overestimate the experimental data. Further, a good agreement was found between the coupled code simulations with the 2D and 3D CFD domain. Figure 7(b) shows the comparison between results obtained adopting both the explicit and implicit coupling scheme, using the same time step (0.005 s) and the 2D CFD domain. The results in terms of LBE mass flow rate and pressure differences between inlet and outlet section of the FPS showed differences lower than 1%. Nevertheless, the use of an implicit scheme allows larger time steps and tends to be more stable than explicit scheme. A sensitivity analysis showed that a time step of 0.025 s (five times greater than the one adopted for the explicit scheme) still allowed us to obtain good results without losing accuracy.Figure 7
Time trend of the LBE mass flow rate for Test 206 in NACIE facility.
(a)
(b)The overall result of the isothermal test is reported in Figure8, in terms of calculated LBE mass flow rate versus the experimental value, averaged over each step of the test. Most of the numerical results stay within an error band of 10% with respect to the experimental results up to a maximum of 12% of discrepancy.Figure 8
Average values of the LBE mass flow rate for Test F400 in NACIE facility.In this application, the main advantage of the coupled calculation was the possibility of calculating local values in the portion of the domain modelled with the CFD code. For instance, the 3D CFD model of the NACIE FPS allowed us to analyse the velocity field inside the test section. Figure9(a) shows the contour plot of the velocity magnitude between the exit of the pin bundle and the outlet section, during Step 6 of Test 206. The region downward the two pins is evidenced by the stagnant and low velocity LBE, caused by the geometrical discontinuity. The latter region is also highlighted in Figure 9(b), where the velocity vectors coloured by axial velocity (along z-axis) are reported.Figure 9
Velocity contour plot [m/s] (a) and vector velocity coloured byz-velocity (b) in the 3D CFD domain, during Step 6 of Test 206.
(a)
(b)


## 3.1. NACIE Facility

NACIE [22] is a Lead Bismuth Eutectic (LBE) loop-type facility. It consists of a rectangular loop made of two vertical stainless steel pipes (2.5′′, Sch. 40) 7.5 m long, acting as riser and downcomer, connected with two horizontal pipes 1 m long. The layout of the facility is shown in Figure 4. A heat source is installed in the bottom part of the riser and consists of two electrical heated rods with nominal thermal power of about 43 kW. A water/LBE “tube in tube” counterflow heat exchanger (HX) is placed on the upper part of the downcomer. The secondary side is fed by water at low pressure. The facility can work in both natural and gas enhanced circulation conditions. Indeed, argon can be injected downstream the heat source, through a 10 mm diameter pipe, housed inside the riser from the expansion vessel.Figure 4
Isometric view of NACIE primary loop.The computational domain for RELAP5/Fluent coupled calculations of the NACIE facility, already described and used in previous work [20], is depicted in Figure 5. It is divided into two nonoverlapping regions: the CFD code was used to simulate the Fuel Pin Simulator (FPS) representing the heat source, whereas the remaining parts of the loop and the HX (both primary and secondary sides) were modelled with the STH code RELAP5. Referring to Figure 5(a),Pipe-130 represents the riser, whileTDV-400, TDJ-405, andBranch-125, at bottom ofPipe-130, simulate the injection of the argon which promote the circulation of the LBE. The expansion vessel (fromPipe-146 toPipe-156) allows the separation between argon and LBE, so that the LBE can flow through the horizontal pipe and then through the HX(Pipe-180) and the downcomer (Pipe-200 andPipe-206). The secondary side water system was modelled by means ofTDV-500,TDJ-505, the HX secondary side annular zone(Annulus-510), andTDV-520.Figure 5
NACIE loop computational domain (a); codes exchanged data (b).
(a)
(b)The FPS was modelled by the Fluent code with both a simplified 2D axial-symmetric domain and a 3D model, shown, respectively, in Figures6(a) and 6(b). The 2D geometrical model was discretized by a structured mesh composed of 7668 rectangular cells uniformly distributed in both the axial and radial coordinates. To model the FPS form losses, mainly due to spacer grids, a constant pressure drop coefficient was assumed in the 2D domain. For this purpose, five distinct interior faces were set as “porous jump,” with each one characterized by an equivalent constant local pressure drop coefficient of 0.7. The 3D geometrical domain was modelled using a symmetry plane, as shown in Figure 6(b). The electrical pins were not included in the domain and heat flux boundary condition (BC) was used at the pin walls. The domain was discretized using 141045 hexahedral elements with refinement close to pins wall, near the inlet and outlet sections. One interior face was set as a porous jump and an equivalent constant coefficient of concentrated pressure drop equal to 0.5 was set in order to simulate the pressure drop due to the spacer grid not simulated in the CFD model.Figure 6
2D (a) and 3D (b) CFD domains of the Fuel Pin Simulator (FPS).
(a)
(b)Concerning the coupled calculation, the exchange of data takes place at the inlet and outlet of the FPS, following the scheme displayed in Figure5(b). Fluent gives the mass flow m˙2 and temperature T2, at the outlet of the FPS, to RELAP5 inTDJ-115 andTDV-112; the pressure p1, obtained from the inlet section of the CFD domain, is set inTDV-110 in the RELAP5 domain. Similarly, the CFD BCs at inlet section, that is, mass flow m˙1 and temperature T1, are obtained from the values computed from RELAP5 in the last volume ofPipe-100, while the outlet surface BC pressure p2 is obtained from the value computed in the first volume ofPipe-120.


## 3.2. Test Matrix

The application to NACIE facility of RELAP5/CFD coupled calculations concerned an isothermal test in gas enhanced circulation (Test 206). During the experimental test, the FPS was switched off and the HX was empty and nonactive. The LBE temperature remained in the range of 200–250°C. The experimental test matrix is shown in Table1, which reports the adopted boundary conditions. Test 206 started with LBE at rest in the facility (0 kg/s mass flow rate); then the gas mass flow controller was switched on to a nominal value of 2 Nl/min (Step 1). The gas flow rate was kept constant for the time needed to have a steady state statistically significant. After that, the other steps came in succession, increasing or decreasing the Ar flow rate with the order delineated in Table 1.Table 1
Boundary conditions of Test 206 in NACIE facility.
T
L
B
E
FPS power
HX power
Nominal Ar flow rate [Nl/min]
200–250°C
0 kW
0 kW
Step 1
2.0
Step 7
8.0
Step 2
4.0
Step 8
6.0
Step 3
5.0
Step 9
5.0
Step 4
6.0
Step 10
4.0
Step 5
8.0
Step 11
2.0
Step 6
10.0


## 3.3. Results

The coupled calculation started with the facility at rest and no gas injected in the riser. After the injection of argon was activated (flow rate to 2 Nl/min), the LBE mass flow rate increased up to a value of about 7.7 kg/s, as shown in the first step of Figure7(a). Increasing the argon flow rate (4-5-6-8-10 Nl/min), the LBE mass flow rate increased up to a maximum of about 14 kg/s. In the second part of the test, gas injection was decreased symmetrically with respect to the increasing ramp and the LBE flow rate followed the same trend. As it can be noticed in Figure 7(a), calculated LBE mass flow rate tended to overestimate the experimental data. Further, a good agreement was found between the coupled code simulations with the 2D and 3D CFD domain. Figure 7(b) shows the comparison between results obtained adopting both the explicit and implicit coupling scheme, using the same time step (0.005 s) and the 2D CFD domain. The results in terms of LBE mass flow rate and pressure differences between inlet and outlet section of the FPS showed differences lower than 1%. Nevertheless, the use of an implicit scheme allows larger time steps and tends to be more stable than explicit scheme. A sensitivity analysis showed that a time step of 0.025 s (five times greater than the one adopted for the explicit scheme) still allowed us to obtain good results without losing accuracy.Figure 7
Time trend of the LBE mass flow rate for Test 206 in NACIE facility.
(a)
(b)The overall result of the isothermal test is reported in Figure8, in terms of calculated LBE mass flow rate versus the experimental value, averaged over each step of the test. Most of the numerical results stay within an error band of 10% with respect to the experimental results up to a maximum of 12% of discrepancy.Figure 8
Average values of the LBE mass flow rate for Test F400 in NACIE facility.In this application, the main advantage of the coupled calculation was the possibility of calculating local values in the portion of the domain modelled with the CFD code. For instance, the 3D CFD model of the NACIE FPS allowed us to analyse the velocity field inside the test section. Figure9(a) shows the contour plot of the velocity magnitude between the exit of the pin bundle and the outlet section, during Step 6 of Test 206. The region downward the two pins is evidenced by the stagnant and low velocity LBE, caused by the geometrical discontinuity. The latter region is also highlighted in Figure 9(b), where the velocity vectors coloured by axial velocity (along z-axis) are reported.Figure 9
Velocity contour plot [m/s] (a) and vector velocity coloured byz-velocity (b) in the 3D CFD domain, during Step 6 of Test 206.
(a)
(b)


## 4. Application to Pool Facility

### 4.1. CIRCE Facility

CIRCE is a liquid metal-cooled pool facility realized and located at ENEA Brasimone Research Centre. It consists of a cylindrical main vessel which can host several kinds of test sections, usually hung from a bolted vessel head. The Integral Circulation Experiment (ICE) test section was designed in order to simulate the primary system of a HLM reactor. It is composed of a Venturi-nozzle flow meter, placed downstream the inlet and connected with the heat source (HS) made of 37 pins arranged in a hexagonal lattice. The upper section of the HS is hydraulically linked to the fitting volume. From this latter, a riser leads the main flow path towards the separator. A nozzle is installed in the lower section of the riser to inject the argon which enhances the circulation of the primary coolant. The separator allows the separation of the Argon from the LBE coming from the riser. The main flow path then goes in the heat exchanger (HX). Another component is the dead volume, placed above the HS, which accommodates the power supply cables. Furthermore, a Decay Heat Removal (DHR) system, hydraulically decoupled from the primary system, is also installed in the pool. A sketch of the facility is displayed in Figure10. More details on the CIRCE-ICE facility and its instrumentation can be found in [23].Figure 10
Layout of the CIRCE facility and ICE test section.The computational domain of the CIRCE-ICE facility was decomposed in two regions: the main pool and the DHR were modelled with ANSYS Fluent CFD code; the ICE test section and the water secondary system of the HX were modelled with the STH code RELAP5/Mod.3.3. Regarding the CFD model of the pool, both a simplified 2D and complete 3D geometrical domains were developed. The 2D geometrical model, shown in Figure11(a), was initially developed with the aim of assessing the methodology employing a limited the computational power. This model assumes the DHR’s axis as axis of symmetry and the flow area in each section is preserved in agreement with the actual facility geometry. The grid is composed of about 242500 elements, mainly structured; the first cell thickness close to the wall is such that y+>30. Nevertheless, the laminar model was adopted for the wall treatment in this work. The 3D CFD domain of the pool and the DHR is shown in Figure 11(b). It represents the ICE test section and the DHR. The grid is composed of 1450000 cells, both with tetrahedral and hexahedral elements. The cover gas above the LBE pool was considered applying the VOF mathematical model in Fluent. RELAP5/Mod.3.3 was adopted for simulating the ICE test section. Figure 11 shows the nodalization of the test section, which models the FPS (Pipe-060 represents the active length), the fitting volume withPipe-090, the riser(Pipe-130), the gas injection system thanks toTDV-003,TDJ-004, andBranch-100, the separator(Pipe-135), and the HX (the primary side isPipe-170). The secondary side of the HX was modelled with three series ofpipes comprised betweenTDV-005 andTDV-006.Figure 11
CIRCE-ICE computational domain for coupled calculations, with 2D (a) and 3D (b) CFD models.
(a)
(b)Concerning the coupling strategy, the application to the pool system CIRCE differs from the application to the loop NACIE. The main difference involves the thermodynamic variables exchanged at the interfaces between the domains. The scheme used for NACIE (Figure5) is unstable when applied to a pool system, since a small variation of the pressure given by RELAP5 to Fluent, at the outlet of the pool (test section inlet), resulted in very large variation of the LBE mass flow rate entering the pool (Fluent domain inlet). A variation in the pressure boundary condition at the pool outlet causes a variation in the LBE level in the pool; even small level difference needs to be adjusted through significant LBE mass flow rate at the pool inlet, mainly due to the large dimension of the pool. Therefore, a new scheme was used according to Figure 11. In this case, Fluent gives to RELAP5 the pressures at both the inlet and outlet sections of ICE, whereas RELAP5 gives the mass flow at the inlet and outlet of the pool (CFD domain). The temperature at the exit of the HX is given to the Fluent pool inlet by RELAP5 and vice versa. Fluent gives the temperature at the inlet of ICE toTDV-010 of RELAP5. This strategy yielded a stable behaviour to the simulation of large pool-type facility like CIRCE, ensuring that the depth from the free surface of the two interfaces is the same for both the domains, so that the pressure due to the column of fluid is correctly considered in both the codes.


### 4.2. Test Matrix

The RELAP5/CFD developed model of the CIRCE facility was employed to simulate the experimental test performed at the ENEA Brasimone and representative of an isothermal test with gas enhanced circulation. The experimental test (Test F400) is composed of several steps characterized by different amount of argon injected inside the riser of ICE test section. During the test, the FPS and the HX were not activated and the LBE temperature was uniform at about 280°C. The boundary conditions of Test F400 are summarized in Table2.Table 2
Boundary conditions of Test F400 in CIRCE facility.
T
L
B
E
FPS power
HX power
Nominal Ar flow rate [Nl/s]
280°C
0 kW
0 kW
Step 1
2.0
Step 7
4.6
Step 2
2.4
Step 8
5.0
Step 3
2.8
Step 9
4.2
Step 4
3.2
Step 10
3.2
Step 5
3.8
Step 11
2.2
Step 6
4.2


### 4.3. Results

In the coupled simulation, the argon was injected in the RELAP5 model through theTDJ-004 (see Figure 11). The gas flow rate was set as boundary condition in agreement with the nominal experimental value. Figures 12 and 13 report the main results of this test. Figure 12 displays the time trend of the LBE mass flow rate for both 2D and 3D coupled calculations. The results are compared with the experimental values and with the RELAP5 standalone calculations. The comparison is quite satisfying in all cases, but the computation with 3D CFD domain better matches the experimental results. In addition, the 3D model provided a more realistic three-dimensional flow pattern inside the pool. Some spikes appeared in the coupled calculations, whereas they did not show up in RELAP5 standalone simulation. The cause of these spikes is related to the choice of the time step; too large time step can cause the loss of accuracy in the numerical results. In this case, the consequence was some oscillations at the beginning of the rapid variation of the gas flow, which takes place in few seconds. Figure 13 shows the average value of the LBE mass flow rate for each step computed by the RELAP5 standalone and by the coupled simulations as function of the experimental values. All the numerical results are set within the error band of 4.0%, as shown in Figure 13.Figure 12
Time trend of the LBE mass flow rate for Test F400 in CIRCE facility.Figure 13
Average values of the LBE mass flow rate for Test F400 in CIRCE facility.The main LBE flow path inside the pool can be examined from the CFD-modelled portion of the domain. Figure14(a) illustrates the contour plot and vectors of the velocity magnitude during Step 8 (argon flow rate 5 Nl/s) in the 2D model. It should be pointed out that the vectors length is normalized and not proportional to its module for graphical reasons. The main flow tends to go downward from the heat exchanger outlet (pool inlet) toward the inlet of the test section (pool outlet). The LBE enters from the heat exchanger outlet (pool inlet), which is modelled through a thin annular section. At the inlet, the fluid spreads in the large volume of the pool, creating a recirculation zone around it. The main flow tends to go down towards the inlet of the test section (pool outlet) at low velocity due to the large flow area. In the lower plenum, the LBE is sucked by the test section inlet. There, the velocity magnitude increases up to 0.8 m/s (full scale range in the zoomed region around the pool outlet), since the flow area reduces considerably. The fluid remains almost stagnant in the upper part of the pool. Some similar considerations can be deduced analysing the CFD results obtained with the 3D domain, shown in Figure 14(b). The LBE enters from the pool inlet, heads for the lower plenum, and spreads in the bottom of the pool. The velocity increases at the pool outlet due to the small cross-section area, achieving a maximum value of about 0.65 m/s. However, the range in Figure 14(b) was set to 0–0.1 m/s to better appreciate the velocity variation in the lower part of the pool. Furthermore, some LBE tends to go from the pool bottom upward and to recirculate in the zone surrounding the FPS pipe. The LBE in the upper part of the pool is almost stagnant. Vectors in Figure 14(b) have an arbitrary length to visualize more clearly the flow direction. Comparing Figures 14(a) and 14(b), it can be noticed that the main flow pattern from the inlet to the outlet is represented in both cases, with comparable velocity magnitude. Nevertheless, some secondary flows as well as recirculation zones are considerably affected by the geometry of the model and just the 3D model is able to predict the main LBE stream realistically.Figure 14
Velocity magnitude field in the 2D (a) and 3D (b) CFD domain of the CIRCE pool, during Step8 of Test F400.
(a)
(b)


## 4.1. CIRCE Facility

CIRCE is a liquid metal-cooled pool facility realized and located at ENEA Brasimone Research Centre. It consists of a cylindrical main vessel which can host several kinds of test sections, usually hung from a bolted vessel head. The Integral Circulation Experiment (ICE) test section was designed in order to simulate the primary system of a HLM reactor. It is composed of a Venturi-nozzle flow meter, placed downstream the inlet and connected with the heat source (HS) made of 37 pins arranged in a hexagonal lattice. The upper section of the HS is hydraulically linked to the fitting volume. From this latter, a riser leads the main flow path towards the separator. A nozzle is installed in the lower section of the riser to inject the argon which enhances the circulation of the primary coolant. The separator allows the separation of the Argon from the LBE coming from the riser. The main flow path then goes in the heat exchanger (HX). Another component is the dead volume, placed above the HS, which accommodates the power supply cables. Furthermore, a Decay Heat Removal (DHR) system, hydraulically decoupled from the primary system, is also installed in the pool. A sketch of the facility is displayed in Figure10. More details on the CIRCE-ICE facility and its instrumentation can be found in [23].Figure 10
Layout of the CIRCE facility and ICE test section.The computational domain of the CIRCE-ICE facility was decomposed in two regions: the main pool and the DHR were modelled with ANSYS Fluent CFD code; the ICE test section and the water secondary system of the HX were modelled with the STH code RELAP5/Mod.3.3. Regarding the CFD model of the pool, both a simplified 2D and complete 3D geometrical domains were developed. The 2D geometrical model, shown in Figure11(a), was initially developed with the aim of assessing the methodology employing a limited the computational power. This model assumes the DHR’s axis as axis of symmetry and the flow area in each section is preserved in agreement with the actual facility geometry. The grid is composed of about 242500 elements, mainly structured; the first cell thickness close to the wall is such that y+>30. Nevertheless, the laminar model was adopted for the wall treatment in this work. The 3D CFD domain of the pool and the DHR is shown in Figure 11(b). It represents the ICE test section and the DHR. The grid is composed of 1450000 cells, both with tetrahedral and hexahedral elements. The cover gas above the LBE pool was considered applying the VOF mathematical model in Fluent. RELAP5/Mod.3.3 was adopted for simulating the ICE test section. Figure 11 shows the nodalization of the test section, which models the FPS (Pipe-060 represents the active length), the fitting volume withPipe-090, the riser(Pipe-130), the gas injection system thanks toTDV-003,TDJ-004, andBranch-100, the separator(Pipe-135), and the HX (the primary side isPipe-170). The secondary side of the HX was modelled with three series ofpipes comprised betweenTDV-005 andTDV-006.Figure 11
CIRCE-ICE computational domain for coupled calculations, with 2D (a) and 3D (b) CFD models.
(a)
(b)Concerning the coupling strategy, the application to the pool system CIRCE differs from the application to the loop NACIE. The main difference involves the thermodynamic variables exchanged at the interfaces between the domains. The scheme used for NACIE (Figure5) is unstable when applied to a pool system, since a small variation of the pressure given by RELAP5 to Fluent, at the outlet of the pool (test section inlet), resulted in very large variation of the LBE mass flow rate entering the pool (Fluent domain inlet). A variation in the pressure boundary condition at the pool outlet causes a variation in the LBE level in the pool; even small level difference needs to be adjusted through significant LBE mass flow rate at the pool inlet, mainly due to the large dimension of the pool. Therefore, a new scheme was used according to Figure 11. In this case, Fluent gives to RELAP5 the pressures at both the inlet and outlet sections of ICE, whereas RELAP5 gives the mass flow at the inlet and outlet of the pool (CFD domain). The temperature at the exit of the HX is given to the Fluent pool inlet by RELAP5 and vice versa. Fluent gives the temperature at the inlet of ICE toTDV-010 of RELAP5. This strategy yielded a stable behaviour to the simulation of large pool-type facility like CIRCE, ensuring that the depth from the free surface of the two interfaces is the same for both the domains, so that the pressure due to the column of fluid is correctly considered in both the codes.


## 4.2. Test Matrix

The RELAP5/CFD developed model of the CIRCE facility was employed to simulate the experimental test performed at the ENEA Brasimone and representative of an isothermal test with gas enhanced circulation. The experimental test (Test F400) is composed of several steps characterized by different amount of argon injected inside the riser of ICE test section. During the test, the FPS and the HX were not activated and the LBE temperature was uniform at about 280°C. The boundary conditions of Test F400 are summarized in Table2.Table 2
Boundary conditions of Test F400 in CIRCE facility.
T
L
B
E
FPS power
HX power
Nominal Ar flow rate [Nl/s]
280°C
0 kW
0 kW
Step 1
2.0
Step 7
4.6
Step 2
2.4
Step 8
5.0
Step 3
2.8
Step 9
4.2
Step 4
3.2
Step 10
3.2
Step 5
3.8
Step 11
2.2
Step 6
4.2


## 4.3. Results

In the coupled simulation, the argon was injected in the RELAP5 model through theTDJ-004 (see Figure 11). The gas flow rate was set as boundary condition in agreement with the nominal experimental value. Figures 12 and 13 report the main results of this test. Figure 12 displays the time trend of the LBE mass flow rate for both 2D and 3D coupled calculations. The results are compared with the experimental values and with the RELAP5 standalone calculations. The comparison is quite satisfying in all cases, but the computation with 3D CFD domain better matches the experimental results. In addition, the 3D model provided a more realistic three-dimensional flow pattern inside the pool. Some spikes appeared in the coupled calculations, whereas they did not show up in RELAP5 standalone simulation. The cause of these spikes is related to the choice of the time step; too large time step can cause the loss of accuracy in the numerical results. In this case, the consequence was some oscillations at the beginning of the rapid variation of the gas flow, which takes place in few seconds. Figure 13 shows the average value of the LBE mass flow rate for each step computed by the RELAP5 standalone and by the coupled simulations as function of the experimental values. All the numerical results are set within the error band of 4.0%, as shown in Figure 13.Figure 12
Time trend of the LBE mass flow rate for Test F400 in CIRCE facility.Figure 13
Average values of the LBE mass flow rate for Test F400 in CIRCE facility.The main LBE flow path inside the pool can be examined from the CFD-modelled portion of the domain. Figure14(a) illustrates the contour plot and vectors of the velocity magnitude during Step 8 (argon flow rate 5 Nl/s) in the 2D model. It should be pointed out that the vectors length is normalized and not proportional to its module for graphical reasons. The main flow tends to go downward from the heat exchanger outlet (pool inlet) toward the inlet of the test section (pool outlet). The LBE enters from the heat exchanger outlet (pool inlet), which is modelled through a thin annular section. At the inlet, the fluid spreads in the large volume of the pool, creating a recirculation zone around it. The main flow tends to go down towards the inlet of the test section (pool outlet) at low velocity due to the large flow area. In the lower plenum, the LBE is sucked by the test section inlet. There, the velocity magnitude increases up to 0.8 m/s (full scale range in the zoomed region around the pool outlet), since the flow area reduces considerably. The fluid remains almost stagnant in the upper part of the pool. Some similar considerations can be deduced analysing the CFD results obtained with the 3D domain, shown in Figure 14(b). The LBE enters from the pool inlet, heads for the lower plenum, and spreads in the bottom of the pool. The velocity increases at the pool outlet due to the small cross-section area, achieving a maximum value of about 0.65 m/s. However, the range in Figure 14(b) was set to 0–0.1 m/s to better appreciate the velocity variation in the lower part of the pool. Furthermore, some LBE tends to go from the pool bottom upward and to recirculate in the zone surrounding the FPS pipe. The LBE in the upper part of the pool is almost stagnant. Vectors in Figure 14(b) have an arbitrary length to visualize more clearly the flow direction. Comparing Figures 14(a) and 14(b), it can be noticed that the main flow pattern from the inlet to the outlet is represented in both cases, with comparable velocity magnitude. Nevertheless, some secondary flows as well as recirculation zones are considerably affected by the geometry of the model and just the 3D model is able to predict the main LBE stream realistically.Figure 14
Velocity magnitude field in the 2D (a) and 3D (b) CFD domain of the CIRCE pool, during Step8 of Test F400.
(a)
(b)


## 5. Conclusions and Perspectives

This work illustrated the main improvements implemented in the coupling methodology, already established at the University of Pisa, between a modified version of RELAP5 code and ANSYS CFD code. A semi-implicit coupling scheme was developed, in which the two codes exchange data within the same time step, until a specified convergence criterion is satisfied. The main advantages introduced by the implicit scheme are represented by an improved stability of the calculation and, consequently, the possibility of adopting larger time steps (compared to the previously developed explicit scheme), which leads to a computational time reduction while preserving the results accuracy. A further improvement consisted in the introduction of Fluent CORBA add-ons to manage more efficiently the Fluent session inside the MATLAB platform through the TUI specific commands.The updated coupling tool was used to simulate a loop system (NACIE facility) and a pool system (CIRCE facility). The application to a pool system required a different set of thermodynamic variables exchanged at the interfaces. Indeed, the methodology used for a loop-type system is unstable when applied to pool-type system. A new stable approach was found and employed in this work.Two experimental tests (Test 206 in NACIE and Test F400 in CIRCE), performed for the hydraulic characterization of the facilities, were reproduced employing the coupling tool to assess the methodology. Results related to the simulation of Test 206 were quite satisfying, even though it tended to overestimate the experimental results with 12% discrepancy maximum. In the author opinion, this discrepancy is linked to inconsistencies in the modelling of the spacer grids, which were not directly modelled but taken into account through local pressure drop coefficients. The comparison of numerical results with experimental data was even more satisfactory in the simulation of Test F400 in CIRCE. In this case, the numerical results approximate the experimental data within an error band 4% of the average experimental value.The main achievement in the use of the developed coupling tool relies on the possibility of exploiting simultaneously the benefits of multiple codes. The required level of detail and/or the need of specific physical models determine the use of CFD or STH codes in different parts of the geometrical domain. The RELAP5 multiphase and multicomponent models are more suitable to investigate two-phase flow regimes (e.g., simulation of the gas-induced flow or the change of phase in the heat exchanger of the secondary water system) and, in general, large portions of the system domain. On the other hand, the CFD code is the appropriate tool to reproduce complex 3D thermohydraulic phenomena. In particular, this work focused on the analysis of the flow path and velocity distribution for both NACIE (pin bundle) and CIRCE (pool) CFD domains, emphasizing the stagnant zones and recirculation areas. The positive assessment reported for this coupling methodology will permit further utilization to more complex nonisothermal conditions. More specifically, the modelling of heat source and/or heat sink will provide, among others, accurate information on temperatures distribution, heat transfer, and thermal stratification phenomena.



---

*Source: 1936894-2017-10-09.xml*