New RELAP 5-3 D Lead and LBE Thermophysical Properties Implementation for Safety Analysis of Gen IV Reactors

The latest versions of RELAP5-3D© code allow the simulation of thermodynamic system, using different type of working fluids, that is, liquid metals, molten salt, diathermic oil, and so forth, thanks to the ATHENA code integration. The RELAP5-3D© water thermophysical properties are largely verified and validated; however there are not so many experiments to generate the liquid metals ones in particular for the Lead and the Lead Bismuth Eutectic. Recently, new and more accurate experimental data are available for liquid metals. The comparison between these state-of-the-art data and the RELAP5-3D© default thermophysical properties shows some discrepancy; therefore a tool for the generation of new properties binary files has been developed. All the available data came from experiments performed at atmospheric pressure. Therefore, to extend the pressure domain below and above this pressure, the tool fits a semiempirical model (soft sphere model with inverse-power-law potential), specific for the liquid metals. New binary files of thermophysical properties, with a detailed mesh grid of point to reduce the code mass error (especially for the Lead), were generated with this tool. Finally, calculations using a simple natural circulation loop were performed to understand the differences between the default and the new properties.


Introduction
The Generation IV International Forum (GIF) is devoted since 2000 to promoting an international collaboration for Generation IV nuclear energy systems development to meet the world's future energy needs.The established requirements for a new generation of nuclear reactors are a more efficient fuel usage, the waste production reduction, the economic competitiveness, and the high standards of safety and proliferation resistance.Following these principles more than one hundred reactor concepts have been examined and six of them have been selected for further research and development.The Lead-Cooled Fast Reactor (LFR) is one of them and it consist in a fast neutron spectrum reactor cooled by liquid Lead Bismuth Eutectic (LBE) or Lead.This concept is particularly interesting because of the good quality of the primary coolant as the low-pressure, chemically inertness and very good thermophysical properties.Because these two fluids are interesting for the next generation of Nuclear Power Plants (NPP) in different safety analysis codes it is just possible to use as working fluid the LBE and the Lead.In the RELAP5-3D © (R5-3D) these two fluids are just available but the relative thermophysical properties do not take into account the latest available experimental data.The Expert Group on Heavy Liquid Metal (HLM) technologies in 2015 published the second edition of the LBE handbook [1] in which new experimental data and updated and/or alternative explanations of phenomena are included.This handbook is now a reference common database in which recommendations relevant to the design of HLM nuclear systems as the Accelerator-Driven System (ADS) and the LFR are described.
In order to align the R5-3D code thermophysical properties for the LBE and Lead to the reference, the Department of Astronautical, Electrical and Energy Engineering (DIAEE) in collaboration with the Idaho National Laboratory (INL) developed a procedure using the software platform MATLAB to generate the new binary file needed by the code.Properties description 4 Description of the fluid (used also in the R5-3D output file), information about the file generation date Triple and critical point information 5 The triple and critical points temperature pressure and specific volume, the admissible pressure, and temperature range for the fluid Number of points' information 6 Information about the number of temperature and pressure points and the total amount of data in block 7 Thermophysical properties data 7 The thermophysical properties data, including the temperature and pressure points This paper provides a description of the procedure and the developed tools used to generate the R5-3D binary file for the implementation of the new thermophysical properties for the LBE and Lead recommended in [1] and is structured as follows.In Section 2 the "tpf" binary file structure is described and the original properties have been compared with the recommended ones, in Section 3 the theory on which are based the soft sphere model and its application to the LBE and Lead is reported, in Section 4 the new "tpf" files properties generation is described, and finally Section 5 compares the natural circulation (NC) system model main parameters using the new and the original properties.

The tpf R5-3D Original Binary File
The thermophysical properties of the fluids used by R5-3D are stored in a number of binary files whose name is composed of the "tpf" prefix and the chemical name of the fluid as suffix (e.g., tpfpb and tpfbipb).Some traces of the binary file structures needed to read and write new binary files are available in [2] for RELAP5/MOD3.3version, but some considerable changes, for example, the implementation of the new XDR format (to avoid machine dependencies of the data), have been implemented in R5-3D version as described in [3].A first MATLAB tool to read the "tpf" files of the RELAP5/MOD3.3was developed during 2014 but this tool is not able to read the R5-3D version of the "tpf" files.
During the last year the efforts to understand the new data structure and the relative updating of the MATLAB tool involved the production of the in-house developed tool "Xfluid" which is able to read and write and postprocess all the data in the R5-3D "tpf" files.

The tpf Thermophysical Properties Binary File Structure.
The "tpf" binary file is divided into seven blocks described in Table 1.At the beginning of each block a keyword and four bytes describe the dimension of the whole block.Using these two headers is possible to parse the binary file and load the data in different MATLAB array for postprocessing purposes.
As described in Table 1, all the thermophysical properties whose distribution is defined in block 6 are included in the last block.The last block is divided into 5 subblocks as described in Table 2.The temperature and pressure grids are included in the first two subblocks, the third and fourth subblocks include 18 vapor and liquid phases properties, calculated in saturation condition along the grids previously defined, finally the last subblock describes 8 subcooled and superheated fluid properties calculated on the points of a bidimensional grid generated using the temperature and pressure points defined in the first and second subblocks.
In Table 2 the superscripts "" and "V" indicate, respectively, the liquid and vapor phase of the fluid, the superscripts "sat" indicate the saturation value of the property, "  " and "  " are, respectively, the total number of temperature and pressure points, "V" is the specific volume (m 3 kg −1 ), "" is the specific internal energy (J kg −1 ), "" is the isobaric thermal expansion coefficient (K −1 ), "" is the isothermal compressibility factor (Pa −1 ), "  " is the specific heat capacity at constant pressure (J kg −1 K −1 ), "" is the specific entropy (J kg −1 K −1 ), "" is thermal conductivity (W m −1 K −1 ), "" is the dynamic viscosity (Pa s), and "  " is the surface tension (J m −2 ).

R5-3D Original Properties versus New Ones Comparison.
The first step of this work was to understand if the original thermophysical properties for the LBE and Lead are different enough from the properties recommended in [1] (in the form of temperature functions), to justify the generation of the new "tpf" files.In [1] there are recommendations only for the liquid phase (for practical application the vapor phase does not appear except for the few saturated vapor in a possible noncondensable gas plenum); therefore a graphical comparison of all the liquid properties was performed using the "Xfluid" tool.In Sections 2.2.1 and 2.2.2 for brevity only the comparison of main properties of the LBE and Lead is reported.The transport properties in the "tpf" files do not show pressure dependencies; instead the thermodynamic properties depend on both the temperature and pressure.For the comparison the data at atmospheric pressure were selected because the correlations recommended in [1] are valid only for this pressure.

The LBE Properties Comparison.
The LBE range of temperatures considered for the comparison is from the LBE melting temperature of 398 K to 1200 K which is a reasonable global interval of validity of the correlation described in [1].This range was selected also because it is adequate for practical application of the LBE properties for safety analysis purposes.The first investigated quantity has been the saturation vapor pressure.The function recommended in [1] came from the experimental data fitting using a modified version of the Clausius-Clapeyron equation.As represented in Figure 1 using the log scale for the -axis, the R5-3D data are several orders of magnitude above the recommended function.Figure 2 shows that the absolute value of the relative error in the selected range of the specific volume is below 2%, but the derivative of the specific volume (i.e., the isobaric thermal expansion) is different.The absolute value of the relative error on the isobaric specific heat in Figure 3 is larger than the 5% near the melting point; it falls below the 5% above the 500 K but remains near this value until 1200 K.In Figure 4 the isothermal compressibility from underestimation of a few % points diverges exceeding the 10%.The thermal conductivity and dynamic viscosity transport properties in Figures 5 and 6 show an average absolute value of the relative error approximately of the 5%, with a maximum, in both cases, near the center of the selected range.Therefore, the relative error is near the acceptability range (1-2%) only for the specific volume.For the other properties the error is >5% except for the saturation pressure, where very large errors have been evidenced and should be reduced.Replacing only the properties with a nonacceptable relative error is not possible because they represent a real fluid in which the physical quantities are correlated; therefore all of them should be recalculated in a consistent way and replaced.Another reason to calculate the new properties is that in the selected range the properties in the "tpf" file have been calculated with a mesh of at least 10 K; therefore it is possible to reduce the mass error during the calculations reducing the mesh as described in [3].

The Lead Properties Comparison. The original Lead
properties were analyzed in the same way.In this case the analyzed range of temperatures was from the Lead melting temperature of 600.2 K until 1500 K, for the same reason described in Section 2.2.1.In this case the saturation pressure data in Figure 7, represented using the log scale for the axis, match very well in a wide range of temperatures with the function recommended in [1], probably because the points in the "tpf" file were calculated using the same function and Relative error (%) the relative error is simply due to the function discretization.Also the specific volume data visible in Figure 8 are in full agreement with the recommended function.In Figure 9 the specific heat at constant pressure match well until 900 K; above this temperature the absolute value of the relative error increases gradually to 4% near 1500 K.The absolute value of the relative error for the isothermal compressibility is higher than 10% in the whole range, as shown in Figure 10.
In Figure 11 the thermal conductivity data show the averaged absolute value of the relative error about 4%, with a maximum above 5% in the center of the selected range.Conversely, as shown in Figure 12, the dynamic viscosity data are in good agreement with the recommended function and the absolute value of the relative error remains below 4%.In this case, several Lead properties are in good agreement with the recommendation, but in the "tpf" file only 52 points for the temperature and 18 for the pressure are available, with a mesh of at least 10 K in the considered range.Therefore also in this case the number of points should be increased, at least to the level of LBE "tpf" file, in which 107 points for the temperature and 105 points for the pressure are included.

The Soft Sphere Model
Unfortunately, an exact Equation of State (EOS) describing the relations between the state variables of a metal in the whole liquid phase existence range is not yet available.Therefore, it is necessary to use some semiempirical models qualified on the available experimental data to extrapolate the thermodynamic properties in the regions where there is a lack of experimental data.One of these models is the modified soft sphere EOS, developed in [4] using the original soft sphere theory described in [5,6] and applying some adjustments for the liquid metals applicability.As reported in [3] the modified soft sphere EOS was already used to calculate the original data in the liquid metals "tpf" files using the available experimental data at the time.Now there is a recommended correlation for each thermodynamic and transport property based on the experimental data available.
The modified soft sphere EOS has been selected again to calculate the thermodynamic properties above and below the atmospheric pressure using the recommended functions to qualify the model.Instead the pressure dependence of transport properties has not been taken into account.This methodology is not a rigorous approach, because in this way the error for the data at different pressures could increase very rapidly, but if it is considered that the liquid metals are almost incompressible fluids and their practical applications are all in the region of the atmospheric pressure, an accurate fit of the recommended functions should ensure the validity of the calculated data.In the following sections a description of the theoretical model and its optimization procedure on the experimental data are presented.

The Theoretical Model Description.
In [4] some attempts were performed to describe the liquid metals behavior in the liquid-vapor and critical regions using the hard-sphere van der Waals model.In this model the fluid particles are defined simply as impenetrable spheres of diameter "" that cannot overlap in space thanks to their pair pure repulsive potential defined as in which "" is the relative distance between two hard spheres.This model has different limitations; therefore its applicability to the liquid metals is not recommended as described in [4].
To improve the theory and experiments agreement a more promising model is the soft sphere EOS.In this model the fluid particles are defined as penetrable spheres which can overlap during collision and the dynamics is determined by the purely repulsive inverse power pair potential law defined in where "" is now the virtual diameter of the spheres, "" is the interaction strength, and "" is the spheres hardness.In [5,6],  Relative error (%)  thanks to an extensive campaign of Monte Carlo calculations using a soft spheres system composed of a total number of particles "" interacting through the pair potential defined in (2), it has been possible to obtain an expression (3) for the potential contribution to the Helmholtz free energy "Φ" relative to a static face centered lattice Φ =    6 (( + 4)  /9 (     ) ) , where "  " is the Boltzmann constant (m 2 kg s −2 K −1 ), "" is the volume occupied by all the particles in the system (m 3 ), "" is the density relative to a closely packed face centered cubic lattice, and "" is the system temperature (K).Equation ( 3) is one of the terms of the partition function "" (5) defined in [6] where "" is Euler's number, "Λ" is the thermal de Broglie wavelength (m), "  " is the chemical specie particle mass (kg), "ℎ" is the Planck constant (J s), "  (/  )" is the attraction contribution term, and "Φ 0 " is the static lattice free energy (J).By using (5) and assuming (7), it is possible to calculate the Helmholtz free energy "" (J) from the partition function through ( 9) ≃ 6.0 + 6.669  − 3 where "  " is the face centered cubic (FCC) lattice Madelung constant, which can be approximated by ( 8) as described in [7].To take into account also the electronic contributions to the heat capacity which may be important for the liquid metal, a nondimensional multiplication factor "" has been added in (3).The last modification to the original model is the introduction of a cohesive energy term " coh " (J) in the final expression of the Helmholtz free energy: where the first term represents the particles kinetic energy, the second term represents the static lattice free energy, the third term is the repulsive potential contribution, and the fourth is the attractive potential energy term.Once the expression for the Helmholtz free energy is available thanks to the well-known correlation between the thermodynamic state variables, it is possible to calculate, for example, the pressure "" (Pa), the enthalpy "" (J), and the isothermal compressibility factor "  " (Pa −1 ) of the system using 3.2.The Optimization Procedure.In (10) the parameters ", , , , " should be adjusted to fit the experimental data available for the selected liquid metal.The standard procedure used in [4,[8][9][10] is based on two steps.First the ", " parameters were calculated by solving system (12) obtained from ( 11) and ( 10) for melting conditions where "ℎ" is the specific enthalpy (J kg −1 ), "" is the molar mass (kg mol −1 ), "  " is the Avogadro number, and " coh " is the specific cohesive energy (J kg −1 ).After the ", " calculation, the remaining parameters ", , " were adjusted to fit the critical parameters and the experimental data.Using the standard approach, the five parameters final values depend on the ", , " guess values and the experimental data fitting accuracy has some limitations because not all the five degrees of freedom were actually used.The new developed procedure uses the MAT-LAB optimization tool.This tool allows finding the minimum of an error function, varying the selected parameters in their domain.To use this tool, it is necessary to define the error function "  (, , , , )" and the parameters constrains which can be linear "4 <  < 15; 0.1 <  < 2, etc." or nonlinear "  (, , , , ) < 0." Generally the error function is defined as where "  " is the experimental data, "(  )" is the model value for the point, "  " is the total number of experimental data, and "  " is the weight associated with the experimental data.
In [1] all the available experimental data have been already analyzed, weighted (basing on the data accuracy), and fitted using temperature polynomials.
To take advantage from the huge work already done on the experimental database, it is possible to define the error function as the integral of the square difference between the temperature polynomials and the model.The specific volume (directly measured), the specific enthalpy (calculated using the measured isobaric specific heat), and the isothermal compressibility (calculated using the data available for the sound speed) are the most important properties which should where " max ,  min " are the bounds of the selected temperature range, "   " are, respectively, the specific volume, specific enthalpy, and the isothermal compressibility temperature polynomials, "  " are, respectively, the model specific volume (obtained solving the first equation of system (12) for atmospheric pressure), the specific enthalpy, and the isothermal compressibility (calculated using (11)), "  " is the property weight based on the importance of the property fitting ( 1 = 100;  2 = 10;  3 = 1), "()" is a weight function which increases the importance of the data near " min ." In the standard procedure after the initial calculation of ", " the others parameters were adjusted to obtain acceptable critical parameters.In Table 3 the critical parameters for LBE and Lead recommended in [1] have been reported with the relative uncertainty.
The critical point is an inflection point for the pressure function; therefore it can be found using system (15) and the first equation of ( 12): Calculating the critical values at each iteration the tool checks if they are in the validity ranges defined in Table 3 or changes the parameters if they are not.Finally, the optimization tool needs some guess value for the five parameters which were calculated imposing the melting condition in (12) and " = 6 ×  0 − 2,  = 1,  = 1" (as recommended in [4]), where " 0 " is the Grüneisen parameter at the normal solid volume.

The LBE Optimization Results.
As shown in Figure 16 using the nonlinear constrain (Con) the critical parameters remain in the validity range (rectangular band) but unfortunately the thermodynamic properties (Figures 13, 14, and 15) do not fit with sufficient accuracy the data and the Normalized Integral Error (NIE) is not acceptable, especially for isothermal compressibility.To obtain better results a new calculation without constrain was performed (noCon).The new calculation largely improves the thermodynamic properties fitting, but only the critical pressure and compressibility are in the validity range.The last calculation was selected to produce the new LBE "tpf" file, because the properties fitting is more important than the critical parameters for the code and in any case the critical parameters assume realistic values.

The Lead Optimization Results
. Also for Lead, as shown in Figure 20, the calculated critical parameters using the constrain are in the validity range.Although this time the critical parameters are exactly on the edge, denoting that the optimization tool does not find a real local minimum, therefore the results of the fitting are very unsatisfactory as shown in Figures 17, 18, and 19.Again, using the optimization without the constrain, the results have been improved significantly.Also in this case the second approach has been preferred for the same reasons described in the previous section and it has been used to generate the new "tpf" file for Lead.

The LBE and Lead New Thermophysical Properties
Once the soft sphere EOS parameters have been optimized for the fluid liquid phase, all the thermophysical properties can be calculated.The work developed in [11] could be used to calculate accurately the LBE vapor phase properties, but it is not relevant for the present calculations; therefore the ideal gas approximation was used for both the LBE and Lead.4.1.Liquid Phase Properties.The thermodynamic properties defined in Section 2.1 have been calculated using the soft sphere EOS and the well-known relations between the properties.In this way the properties pressure dependence is conserved and the specific entropy, the isobaric thermal expansion, and so forth, which were calculated integrating 600 700 800 900 1000 1100 1200 1300 1400 1500 0.5 or differentiating other thermodynamic properties, are completely coherent.The transport properties values have been calculated directly using the functions recommended in [1].

Vapor Phase Properties.
The vapor phase thermodynamic properties as described in [12] have been calculated  using the simple ideal gas model using a constant value for the isobaric specific heat and the well-known correlation for the other properties.The conductivity and the viscosity have been calculated using ( 16) and (17), respectively, as described in [13] where "V  " is the solid molar volume (m 3 mol −1 ) at the melting temperature "  " (K).Finally the integral quantities as the specific internal energy and the specific entropy have been calculated using as starting point the saturated liquid respective value plus the latent heat of boiling.

The Natural Circulation System Model.
To test the new generated properties, a simple NC system has been nodalized.The main parameters of the model are reported in Table 4.These parameters have been selected to maintain the fluid temperature at the steady state and during the transient between 650 K and the 800 K for both liquid metals.As shown in Figure 21 the heater is located at the bottom of the hot leg and the heat exchanger is located at the top of the cold leg to maximize the driving head.The heat exchanger consists in a simple "AISI-316" plate with the NC system on the left side and a fixed wall temperature on the right side.The expansion tank is connected to the highest point of the loop to allow the fluid expansion through an Argon cover gas plenum.The liquid metal level is initialized at the middle height of the tank and the pressure is imposed from the top with a time dependent volume.A 90 ∘ elbow pressure drop coefficient is introduced at the boundaries of the hot leg and at the end of the cold leg; a "" junction pressure drop coefficient has been introduced at the beginning of the cold leg to take into account the expansion tank connection.

Results
Comparison.To test the numerical stability of the original and new properties, a first steady state simulation of 800 seconds with a time step of 1 ⋅ 10 −3 s has been performed.The system estimate mass error calculated using the new properties (Figure 22) is lower for both liquid metals.Also the truncation mass error fraction (Figure 23) calculated as described in [3] using the new properties appears to be lower especially for the Lead.These achievements derived from the increasing of the number of temperature and pressure points in the new "tpf" files.
After the steady state calculations, a transient analysis has been performed to test the new properties in different operative ranges.The transient consists in three power ramps each one of 80 kW in 50 s, followed by 800 s of steady state to   stabilize the system.In the following sections some significant results are presented.

LBE Results
Comparison.One of the most important parameters in a NC loop is the mass flow rate, whose value is determined by many different fluid properties.Therefore it should be a good parameter to estimate the effect of the differences between the original and the new thermophysical properties.As shown in Figure 24, the mass flow rate increases in the calculation with the new properties due to the increase in the temperature (and density) difference between the two legs, thus leading to a higher mass flow rate.The heat transfer coefficient as shown in Figure 26 is underestimated by about 5%, but the error shows a downward trend during the transient.

Lead Results
Comparison.The original Lead thermophysical properties overestimate the mass flow rate and the system irreversible pressure losses, as shown in Figures 27 and 28, but the error is quite limited, about 0.5% and 0.7%, respectively.The heat transfer coefficient, shown in Figure 29, starts instead from an overestimation of about 0.4%, rising up to 2% at the end of the transient.

Conclusions
In this paper a full comparison between the original thermophysical properties of R5-3D and the NEA recommended ones has been reported for two heavy liquid metals, LBE and Lead.To tests the discrepancy effects on the calculations, an in-house developed tool "Xfluid" has been developed to generate new "tpf" files from NEA recommendations using the soft sphere model, whose theoretical model and the adopted optimization procedure on the experimental data have been described in detail.A simple NC system model has been used to perform some comparison calculations.The NC system main parameters calculated using the new and original LBE properties show a significant discrepancy; therefore further investigations and validation using experimental data should be performed.On the contrary the parameters calculated using the new and original Lead properties show a limited discrepancy, except for the heat exchange coefficient, whose error increases during the transient.Furthermore, even if the main parameters are globally in good agreement for Lead, the truncation mass error using the original "tpf" file is higher than the other cases: therefore it has been reduced increasing the temperature and pressure points in the new Lead "tpf" file.The future activities will be devoted to find an optimal pressure and temperature grid to minimize the numerical issues during the calculations and to the validation of the new thermophysical properties using experimental data.

Figure 18 :
Figure 18: Lead specific enthalpy, NEA data versus constrained and nonconstrained optimization results.

Figure 22 :
Figure 22: System mass error, LBE and Lead original versus new properties.

Figure 23 :
Figure 23: System truncation error, LBE and Lead original versus new properties.

Table 3 :
The LBE and Lead critical parameters.

Table 4 :
Natural circulation system main parameters.