Preliminary CFD Assessment of an Experimental Test Facility Operating with Heavy Liquid Metals

The CFD analysis of a Venturi nozzle operating in LBE (key component of the CIRCE facility, owned by ENEA) is presented in this paper. CIRCE is a facility developed to investigate in detail the fluid-dynamic behavior of ADS and/or LFR reactor plants. The initial CFD simulations have been developed hand in hand with the comparison with experimental data: the test results were used to confirm the reliability of the CFD model, which, in turn, was used to improve the interpretation of the experimental data. The Venturi nozzle is modeled with a 3D CFD code (STAR-CCM+). Later on, the CFDmodel has been used to assess the performance of the component in conditions different from the ones tested in CIRCE: the performance of the Venturi is presented, in terms of pressure drops, for various operating conditions. Finally, the CFD analysis has been focused on the evaluation of the effects of the injection of an inert gas in the flow of the liquid coolant on the performance of the Venturi nozzle.


Introduction
Among the GEN-IV proposed concepts [1], the gas cooled [2,3] and Heavy Metal (HM) cooled [4] reactors seem to be particularly interesting, also for fuel cycle closure purposes [5,6].Regarding HM reactors, many experiments have been and are carried out in order to improve the knowledge of the coolant behavior and to validate the simulations tools (codes).CIRCE (Circulation Experiment, shown in Figure 1(a)) is a large-scale test facility designed for studying key operating principles of the 80 MW Experimental Accelerator-Driven System (XADS) located in Brasimone, Italy.The activities planned for this facility aims at investigating the hydraulic, chemical, and mechanical issue related to the use of the eutectic mixture of lead and bismuth (LBE) as a coolant, in a pool configuration.The test facility was extensively instrumented and tested in order to reliably perform the subsequent campaign of thermal-hydraulic experiments [7].
CIRCE plant is filled with LBE, which is a eutectic alloy of lead (44.5%) and bismuth (55.5%) used as a coolant in some nuclear reactor and in this facility, too.One of the key features of LBE is that the melting temperature (125 ∘ C) is definitely lower than the one of pure Pb (327 ∘ C).This fact is particularly useful because it reduces heating requirements and the risk of unintended solidification of the coolant.Moreover, at lower temperature, the corrosion rates of materials are significantly reduced and so this allows the use of conventional code-qualified materials for very long test durations.A disadvantage is the neutron capture of Bi 209 and its subsequent beta decay that form Po 210 which is an  emitter.One of the main drawbacks of using LBE (or, similarly, pure lead) as coolant is that it has a strong effect of corrosion and erosion on steel, which limits the range of the velocity and of the temperature of the coolant [8].
In this frame, it is particularly interesting to evaluate the performance of the static devices, such as jet pumps (Figure 1 and low-Reynolds numbers.Jet pumps and ejectors are static pumping devices, which exploit the pressure and momentum of a little mass flow rate to move, and pump, a larger one.The effectiveness of the jet pump directly depends on the process of mixing of the two confined jets, which regulates the way in which the fast jet transfers its momentum and its kinetic energy to the entrained slow fluid [9].

Flow Equations and Numerical Model
The classical form of Reynolds-Averaged Navier-Stokes equations stands for liquid metals, too, as far as their rheological behavior corresponds to the Newtonian fluid model.The governing equations of fluid flow represent mathematical statement of the conservation laws of mass and momentum for isothermal flow, and, in the case of stationary motion, the form is as follows: Mass conservation is Momentum equation is ( The thermal conditions of the test facility were carefully controlled to prevent the solidification of the mixture; therefore, in the Venturi nozzle, which is placed inside of the facility, heat transfer to the surroundings is negligible and the temperature is constant through the fluid domain.For this reason, the energy equation has not been solved and the flow is assumed isothermal and adiabatic.
The turbulence is treated by a RANS approach.The RANS models implemented in STAR-CCM+ [10] and used for these simulations are as follows: The turbulent boundary layer can be divided into three zones: viscous sublayer, log layer, and outer layer [11].
STAR-CCM+ uses a wall treatment to deal with the turbulent boundary layer.As is known, a wall treatment is the set of near-wall modeling assumptions for each turbulence model.
Three types of near-wall treatment are provided [10]: (i) The high- + wall treatment implies the wall-functiontype approach.It is suitable for use with models that do not explicitly damp the turbulence in the near-wall region.
(ii) The low- + wall treatment is suitable for low-Reynolds number turbulence models.It is possible to use this treatment only if the entire mesh is fine enough.
(iii) The all- + wall treatment is a hybrid treatment that attempts to emulate the high- + wall treatment for coarse meshes and the low- + wall treatment for fine meshes.It is also formulated with the desirable characteristic of producing reasonable answers for meshes of intermediate resolution.

CFD Analysis of the Venturi Nozzle
One of the studies carried out at Brasimone Research Center is the calibration of a large Venturi nozzle flow meter (Figure 2) operating in LBE [12].Such flow meter has been selected to provide flow rate measurements during thermalhydraulic test on eutectic Pb-Bi (LBE) fluid.Because of the features of this fluid, and the importance of the tests, the Venturi nozzle underwent to a specific calibration aimed at the definition of the characteristic response of the device, with LBE as working fluid [13].The aim of this calibration is to define the correlation among the differential pressure evaluated between the inlet and the throat, and the Pb-Bi mass flow evolving in it.
As already anticipated, the CFD software used in these simulations is STAR-CCM+ by CD-Adapco [10].The 3D-CAD tool provided by this software allows modeling the fluid domain.As shown in Figure 3, the geometry was created with a 360 ∘ revolution of a 2D sketch.
The boundary conditions set for the simulations are as follows: (i) Mass flow inlet: imposed mass flow rate [kg/s] (ii) Pressure outlet: 0 [Pa rel ] (imposed) (iii) Wall: no-slip condition In any flow condition of the present paper, the Reynolds number is always greater than 4000 (the lower limit for turbulence regime), so a turbulent approach is fully justified.
The first set of simulations were devoted to a sensitivity analysis of the effect the mesh, of the turbulence models, and of the variation of the physical properties caused by different operating temperatures.
Figure 4 plots the results obtained with various mesh sizes and realizable - turbulence model: this figure reports in abscissa the ratio between the calculated pressure drop and the dynamic pressure at the inlet of the throat section and in ordinate, Figure 4(a) the ratio between the design flow rate,  nom (equal to 2.25 − 02 m 3 /s), and two different off-design conditions, due also to density variations of the operating fluid, and Figure 4(b) the grid size in kNodes (logarithmic scale).The starting point for the base size (grid dimensions) was a value of 0.02 m, halving this parameter for each simulation reaching the final (and much heavier in terms of computational effort) value of 0.0025 m.
The sensitivity analysis of the mesh size shows that results are fairly grid independent; that is, that the solution would not be significantly influenced by a further mesh refinement.
It was found that the mesh with about 70000 cells and a mesh size of 0.01 m gives the best compromise between calculation time and precision.The wall treatment used is the All- + provided by STAR-CCM+ and the  + values are between 1.5 and 3.1.
Under the hypotheses of steady, inviscid, incompressible fluid, adiabatic flow and negligible change in height, the Bernoulli's equation (see ( 3)) was used for the estimation of the pressure drop among the throat: The results of this simple evaluation are in very good agreement with the experimental data (Figure 4).
As a second step, different turbulence models have been considered.The models used were a realizable -, - standard, and SST.The results obtained by varying the turbulence models are substantially the same obtained with - models, as Figure 5 shows (percentage difference lower than 20%).After these preliminary analyses, a series of new tests have been performed in order to consider the variation of Venturi nozzle correlation with different nominal temperature (200 ∘ C, 300 ∘ C, and 385 ∘ C) and consequently different densities, dynamic viscosities, and mass flow rates at the inlet.For each temperature, three values of flow rates (a minimum value equal to 0.015 m 3 /s, a nominal value equal to 0.0225 m 3 /s, and a maximum value equal to 0.035 m 3 /s) were chosen to compare the different tests.As Figure 6 shows, the pressure drop and mass flow rate values are quite similar so it can be concluded that the effect of temperature is negligible, at least in the investigated temperature range: so this justifies the adoption of a single best fit line.
According to the obtained results, Figure 7 shows the trend obtained with the STAR-CCM+ simulations.This trend is compared with experimental data previously obtained by ENEA [12] and with the inviscid results obtained by Bernoulli's equation.It can be noted that the simulations with STAR-CCM+ provide valid results comparable with the experimental data although a certain deviation exists,  especially with high value of flow rates; the parabolic trend, compared with the Bernoulli's results, confirms the validity of the data.The reason of the almost constant deviations in the range of mass flow rate tested could be due to several factors: in particular, real roughness value of the experimental apparatus, material properties adopted and their comparison with the real properties of the liquid metal flowing in the test facility and, looking at numerical issues, the relative low experience with numerical solver for liquid metal fluid.

Analysis of Different Geometries.
The preliminary CFD analysis of the Venturi nozzle operating condition corresponding to the experimental ones led to the final selection of a reliable mesh.On this basis, the next step has been the evaluation of possible improvements of the nozzle shape, which could lead to better performance as well as to simplification of the device design.Particularly the analysis has been focused on the convergent and on the divergent parts of the nozzle.The implemented modifications are summarized in Figure 8 and include the following: (i) The replacement of the fillet in the inlet throat section by a 45 ∘ chamfered corner (ii) The variation of the diffuser angle and therefore of the diffuser section length Firstly, we dealt with the effect of modifying the convergent part of the nozzle.

Fillet versus Chamfer.
Further analysis concerning different geometry for the nozzle was conducted in order to assess how the geometry influences the flow inside the nozzle.The original inlet throat is characterized by a fillet, where radius is equal to 3.26⋅10 −2 m, to ensure a smooth distribution of the fluid in that section and to reduce flow recirculation and detachment due to section variations and insufficient turbulent mixing, due to the large fluctuation of velocity distribution.The modified geometry analyzed in this case presents a 45 ∘ chamfer instead of the fillet.The results are presented in Figure 9.The different distributions of velocity  inside the Venturi for the two geometries are also presented (Figures 10 and 11).In Figure 11, the jet-like behavior in the throat can be seen, with a higher velocity value and a wider recirculation zone and displacement thickness near the wall region of the straight section next to the throat, with the consequent smaller effective cross section.Due to the nature of the considered Venturi (i.e., a measuring instrument), the region of higher turbulence with creation of vortex shedding nearby the measurement region with the chamfered geometry could lead to inconsistent or not very accurate measures; another problem due to the vortex shedding is the high nonlinearity of the characteristic trend.The amplified pressure drop across the throat (virtually better in terms of accuracy of the results) is not convenient enough due to the previously exposed reasons.
Figure 12: Pressure drop across the diffuser.

Diffuser Analysis.
Regarding the diffuser, we tried to change its total length and vary the angle to observe how these parameters influence the velocity field through the outlet tube, preserving their total length for installation constraints.
Defining the section of the Venturi meter as in Figure 8, it is possible to define [14] a diffuser pressure recovery coefficient and, as mentioned in [14], for uniform flow at the entrance and at the exit and neglecting frictional losses, the ideal static pressure recovery coefficient is defined, for this case, as and, finally, the theoretical efficiency as As shown in Figure 12, where in the abscissa the values obtained by (4) are reported, the diffuser geometry design provided by constructor has the best value by the point of view of pressure losses and a quasi-linear trend in the analyzed range.Low pressure drop across the Venturi was obtained with a high-pressure recovery through the diffuser and a low friction and entropy (due to the vortex shedding development) pressure loss.
Figures 13, 14, and 15 highlight the differences among the considered diffuser geometries.

Transient Simulations.
The Venturi nozzle is one of the components of the control chain of the CIRCE facility.For this reason, it is useful to analyze the dynamic behavior of the nozzle in mass flow rate transient condition.For this reason, the final set of simulations is related to a transient regime: the aim of this activity is to predict the response of the Venturi when mass flow rate, at inlet, varies with time.The trend assumed as inlet boundary condition is taken from [12].
The simulation period has a duration of 100 seconds because it is assumed that smoother variations do not lead to unsteady behavior of the Venturi nozzle.Due to the nature of the flow motion driving force, that is, to the hydrostatic pressure difference caused by free-level differences, the trend is semilinear only in the 100 seconds analyzed; in the final part, the pressure drops down very quickly, due to the level equalization and to the very small fluid velocity.
For these simulations, the density and viscosity values have been calculated for a temperature of 300 ∘ C. The time discretization adopted for this analysis is an implicit method with a time step of 0.01 seconds and the number of iterations  In the first seconds of the simulation (not reported in Figure 16), there is a fluctuation due to an adjustment of mass flow rate and pressure.Then the values tend to decrease according to mass flow rate law imposed as the Venturi inlet condition.The results obtained with STAR-CCM+ are compared with experimental data [12] in the Figure 16 (in the same figure the results obtained by ENEA with FLUENT [12] are also reported).
It is quite evident the difference between the experimental and both numerical results obtained with ANSYS FLUENT and STAR-CCM+: in the first case, the pressure drop is overestimated conversely to the STAR-CCM+ underestimation with respect to the experimental results.The reasons of the underestimation could be assumed due to the lack of a detailed knowledge of effective roughness and LBE physical proprieties during the experimental setup, as well as the linearization of the mass flow trend, which affects the possible generation of entropy due to chaotic behavior of the flow with accumulations and acceleration/deceleration local trends.

Multiphase Analysis
The last set of simulation was aimed at defining the fluiddynamic behavior of nozzle in case of two-phase flow.
As known, reactivity measures the degree of change in neutron multiplication in a reactor core and is related to the tendency of the reactor core to change power level.The reactivity void coefficient is a parameter that can be used to estimate how much the reactivity of a nuclear reactor changes as voids (typically gas bubbles) are growing in the reactor.
In a LFR, boiling is not considered as a relevant scenario due to high boiling temperatures of lead (about 1743 ∘ C) or LBE (about 1670 ∘ C) well above the temperature level that implies that the integrity of the reactor is lost.
Instead, two-phase flow may come from presence of argon in the reactor [15].If the gas diffuses in the coolant and reaches the core, or, due to the presence of a free surface between argon and lead in the pool of the reactor, if the free level is not stationary, for reason like seismic oscillation or reason related to the design of the primary system, then there can be waves on the free surface and a small quantity of argon could be carried by the lead and, finally, can flow from the pump to the core.Argon in the core could cause a positive reactivity insertion.Therefore, it is very interesting to evaluate the possibility of detecting the presence of argon in the coolant by means of measurements of the changes in pressure drops through the Venturi nozzle.
In this analysis, we simulated a multiphase flow of LBE and argon.The gas fills the volume between LBE interface and roof of a LFR, as is further explained in [16].
STAR-CCM+ provides different models to meet the requirements of multiphase flows; in this case, the volume of fluid (also called VOF) has been used [10].This model is provided for systems containing two or more immiscible fluid phases, where each phase constitutes a large structure within the system [17].
In the multiphase simulations, the same model for geometry and mesh used in Section 3.1 has been adopted.The two phases (i.e., argon and LBE) were defined separately with their own properties.An additional parameter, the volume fraction of each phase, defined as the volume of each fluid (LBE or argon in this case) divided by the total volume of the mixture with the constraint indicated in (8), is necessary to define initial fluid characteristics.For this preliminary analysis, assuming very conservative values in order to initially estimate the Ar trapping, this parameter has been varied from a minimum value of argon equal to 0.02 (2%) to a maximum value equal to 0.1 (10%).In the volume of fluid approach, the equations are solved for an equivalent fluid whose physical properties are calculated as functions of the physical properties of its constituent phases and their respective volume fraction.
(i) Conservation of mass for phase  [11] is (ii) Constraint on volume phase is (iii) Material properties (density in the example) is (iv) Single momentum equation is solved throughout the domain, and the resulting velocity field is shared among the phases; these equations, shown below,  are dependent on the volume fractions of all phases through the properties  and : where (i)  is the volume fraction; (ii)  and  is, respectively, the density and the dynamic viscosity of the fluid or the mixture; (iii) the subscript ,  was referred to at the two different fluids or phases; (iv) ṁ is the mass transfer from phase  to phase  and conversely; (v)   . is the mass source terms, equal to zero by default; (vi) V is the velocity of the phases or of the mixture; (vii) the terms  and  in (10) are the gravitational and the external body forces.
The energy equation is not considered because there is no heat exchange with external and the phases have the same temperature.In this steady-state analysis, we also choose to perform simulations with different values of mass flow rate in order to observe how the velocity and the pressure change in the Venturi nozzle.The results obtained are reported in (Figures 17-21).The phenomenon is very complicated, because there are many parameters to be considered.
(i) The insertion of gas phase in a liquid phase causes a local density decrease resulting in the establishment of a upward thrust, with a consequent pressure recovery contribution: assuming that the density in the inlet section can be   evaluated by the homogeneous model, thereby introducing the volume fraction, and introducing the buoyancy force, It is possible to evaluate the "virtual" differential pressure: where ℎ is the height of the section.With this simple analytical estimation of the order of magnitude of the upward trust contribution due to the density reduction, it was noted that it is comparable with the differential pressure due to the Venturi geometry.For example, assuming the values reported in Table 1 and considering in the nominal case a calculated pressure drop about 4000 Pa, it is clear that the buoyancy effect is not negligible.
(ii) There is a different behavior of the pressure drop losses and recovery in the convergent and divergent part of the Venturi, with respect to the single-phase flow.
In the following, only a phenomenological analysis of the phenomenon has been performed with the aim of investing the applicability of the Venturi meter like an argon detector.
Figures 17-19 show the pressure profiles, at the centerline of the Venturi meter, with different value of the Ar volume fraction: it is possible to remark how the pressure drop in the convergent nozzle was approximatively the same for all cases, and the difference could be comparable with the uncertainty of the available (existing) meter.The buoyancy phenomenon is very amplified at low mass flow rate, while at high mass flow rate is less noticeable.
In Figure 20, the pressure drop/recovery across the Venturi nozzle is reported, for different value of mass flow rate and Ar volume fraction: for value larger than 0.05 in Ar volume fraction, the Venturi works like a pump, due to the previously exposed buoyancy phenomenon.In particular, with the lowest mass flow rate, a small insertion of Ar is enough for a change in the pressure drop sign.
The percentage difference between the single-phase and the multiphase differential pressure is reported in Table 2.The reference value for evaluating the relative difference is the differential pressure measured in the single-phase analysis.It is possible to note that, in particular for low and high flow rates, the differential pressure measured varies with the increasing in gas volume fraction.

Transient Analysis.
A final evaluation to complete the multiphase analysis has been carried out in a transient regime.These simulations have been carried out in the same way as explained in Section 3.2.Figure 21 shows the results.As it can be noted, the pressure drop across the Venturi nozzle seems to be a good parameter, also in transient cases, to evaluate the presence of argon in LBE.As it is shown in the steady state, an increase in the value of Ar volume fraction causes a greater decrease of pressure drop across the Venturi nozzle.

Conclusions and Future Developments
The results obtained from the CFD analysis of the Venturi nozzle flow meter of CIRCE facility bring us to these conclusions: (i) The 3D CFD model adopted has provided valid results comparable with the experimental data although a certain deviation is present.
(ii) The comparison between different geometries for the throat inlet shows why the fillet solution should be preferred to the chamfer one.
(iii) The STAR-CCM+ simulations for the unsteady analysis show good results (if compared to experimental trend) for pressure drops.
(iv) Looking at multiphase analysis, a decrease in pressure through Venturi nozzle is observed, depending on argon volume fraction, with respect to single-phase Science and Technology of Nuclear Installations flow.This behavior has been confirmed also by the transient results.
(v) An important remark is the capability of this flow meter to detect the presence of Ar in the lead flow by monitoring the total Venturi pressure drop instead of the throat pressure drop.
In terms of future developments, we can preliminarily say the following: (i) In this work, the roughness was not considered properly.
(ii) Regarding transient simulations, STAR-CCM+ underestimates the experimental trend for pressure drops, while FLUENT overestimates them.It will be also interesting to investigate deeper this difference between the two codes.
(iii) Regarding multiphase, it will be interesting to investigate and (possibly) find a correlation between argon volume fraction and pressure drop.This could allow the use of Venturi nozzle to detect the presence of argon in the flow by simply measuring the pressure drop.
(iv) A phenomenological investigation of the buoyancy contributions in the pressure drop across the Venturi meter will have to be carried out, in order to define the parameters (e.g., lead velocity, Reynolds number, and Venturi geometry) that most affect the considered phenomena.
(v) Finally, although a homogeneous model has been adopted, it could be interesting to further investigate the reasons behind the apparently anomalous pressure increase near the inlet at higher void and its relation with buoyancy over pressure drop ratio.

Figure 1 :
Figure 1: (a) CIRCE facility plant layout and (b) schematic representation of the tank with the indication of the Venturi test section [12].

Figure 5 :
Figure 5: Comparison between different turbulence models and experimental results.

Figure 8 :
Figure 8: Modifications of the Venturi meter geometry.

Figure 9 :Figure 10 :Figure 11 :
Figure 9: Comparison between fillet and chamfer solutions for the inlet throat.
, Cp Design di user, Cp Long di user, Cp Short di user,  Design di user,  Long di user,  Cp,  [-]

Figure 16 :
Figure 16: Comparison between pressure drops obtained from experimental data, STAR-CCM+, and FLUENT simulations.

Figure 17 :
Figure 17: Pressure distribution at minimum mass flow rate for different Ar volume fractions.

Figure 18 :
Figure 18: Pressure distribution at nominal mass flow rate for different Ar volume fractions.

Figure 19 :
Figure 19: Pressure distribution at maximum mass flow rate for different Ar volume fractions.

Figure 20 :Figure 21 :
Figure 20: Pressure drop across Venturi meter for different mass flow rate and Ar volume fractions.

Table 1 :
Pressure difference due to the buoyancy estimation.

Table 2 :
Percentage differences for multiphase analysis.