Nuclear Engineering and Design 359 (2020) 110482
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
Numerical and experimental campaigns for lead solidification modelling and testing
T
Manuela Profira, , Vincent Moreaua, Tomáš Melicharb ⁎
a b
CRS4, Scientific and Technological Park, Ed. 1, 09010, Pula, Cagliari, Italy Research Centre Rez, Hlavni 130, Husinec-Rez 250 68, Czech Republic
ARTICLE INFO
ABSTRACT
Keywords: CFD HLM Pool modelling Thermal hydraulics Lead solidification
The Computational Fluid-Dynamics (CFD) modelling of Heavy Liquid Metal (HLM) flows in pool configuration is investigated in the framework of the SESAME project. This paper focuses on the coupling between numerical simulations and experimental activity with the objective to make CFD a valid tool in support to the design of safe and innovative Gen-IV nuclear reactors. The attention is focused on the possible occurrence of lead solidification phenomenon. The aim is to demonstrate that the overall modelling of HLM complex systems can include solidification phenomenology by reproducing small/medium scale experiment, making comparison with experimental results, improving the numerical setting (post-test) and evaluating time scales, limitations and computational costs. A dedicated experimental campaign on the SESAME-Stand facility has been performed by the Research Centre Rez (CVR), relying on the use of CFD support, with the specific objective to build a series of datasets suited also for the CFD modelling validation. The numerical simulations of the SESAME-Stand experimental facility performed in STAR-CCM + are described. The model uses liquid lead as working fluid in the pool and air in the cooling channel. By increasing the mass flow rate in the cooling air channel, the solidification process is initiated and the freezing front propagates in the pool until reaching an internal obstacle. The issues encountered in the pre-test simulations have been fully overcome by means of systematic investigations and all the improvements have been applied in the post-test model. A decisive modelling aspect was the correct implementation of the thermal radiation which plays an important role in the cooling process. The numerical model allowed to reach an increased number of initial steady states and accessible transients, according to the experimental matrix. The comparison with the experimental results shows similar temperature configurations and comparable lead frozen fractions in the steady state cases and good agreement in the fast transient solidification-remelting cases.
1. Introduction The development and the validation of computational models able to simulate the liquid metal thermal-hydraulics and the relative freezing phenomena are essential for the development of liquid metal cooled Gen IV nuclear reactors and for their safety assessment. The coolant solidification is considered one of the basic phenomena in the context of the thermal-hydraulics challenges (Roelofs et al., 2016). The coolant solidification in liquid metal cooled reactors is a potential safety issue, since it may occur during shut-down, maintenance, refueling and whenever an overcooling from the secondary system or the (passive)
emergency cooling system is verified. However, solidification might be expected even in normal operation in case of occurrence of flow stagnation zones near the cooled surfaces. The local solidification of the liquid metal will affect the natural convection pattern in the pool leading consequently to core overheating. Solidification models are often used in the automotive or metal molding industries. However, interaction of solidification interface with local flow conditions is especially important for the study of liquid metal nuclear reactors. Nevertheless, applicability of such models to the nuclear reactor context has not yet been investigated. Liquid metals considered in reactor applications have melting temperatures above the
Abbreviations: ALFRED, Advanced Lead Fast Cooled European Demonstrator; CAD, Computer Aided Design; CFD, Computational Fluid Dynamics; CFL, CourantFrierichs-Lewy; CIRCE, Circulation Eutettico; CRS4, Center for Advanced Studies, Research and Development in Sardinia; CVR, Research centre Rez; HS, Heat Source; HF, Heat Flux; HT, Heat Transfer; HTC, Heat Transfer Coefficient; MFR, Mass Flow Rate; MYRRHA, Multi-purpose hYbrid Research Reactor for High-tech Applications; NRG, Nuclear Research and Consultancy Group; SESAME, Thermal hydraulics Simulations and Experiments for the Safety Assessment of Metal cooled reactors; TC, Thermocouple; VoF, Volume of Fluid ⁎ Corresponding author. https://doi.org/10.1016/j.nucengdes.2019.110482 Received 20 June 2019; Received in revised form 12 December 2019; Accepted 12 December 2019 0029-5493/ © 2019 Elsevier B.V. All rights reserved.
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Nomenclature φ σ e T Tbottom Tair-In Ttop Tair-Out TEW α T1 T0
T l k0 k+ h q Cp r μ κ ρ
Heat flux [W/m2] Stefan-Boltzman constant [W/m2-K4] emissivity Temperature [C or K] Vessel bottom temperature Air channel bottom temperature Vessel top temperature Air channel top temperature Air channel external wall temperature View factor Top argon layer temperature Bottom argon layer temperature
t
ambient temperature. In certain accidental scenarios with overcooling in primary heat exchangers or through the vessel walls, the liquid-metal coolants can solidify. This might result in partial or complete blockage altering the coolant flow paths in the reactor. Solidification might also lead to mechanical stresses in components as a result of thermal contraction and expansion. Hence, it is important to know if and where solidification can occur during accidental transients and how the solidification front moves. In the current Liquid Metal Fast Reactors (LMFR) framework, activities are undertaken to expand the knowledge on solidification and its propagation in typical liquid-metal reactor coolants (Roelofs et al., 2016). Fundamental experiments are performed, providing validation data for numerical models of solidification. Numerical models need to be further developed and/or validated in parallel. For this purpose, activities focused on development and validation of computational models simulating the liquid lead natural convection and freezing have been performed within the Horizon 2020 SESAME project (SESAME, 2019). The larger context of the pool modelling includes numerical works dedicated to various experimental facilities. After the extended experience gained with the modelling of the entire primary loop of MYRRHA (Moreau, 2014); (Koloszar and Moreau, 2016); CRS4 and NRG dedicated the attention to the modelling of CIRCE (Zwijsen, et al., 2018); (Zwijsen et al., 2019), each partner dealing with complementary tasks, according to a redundancy and diversification criteria. The acquired know-how was extrapolated to the modelling of the ALFRED demonstrator (Visser et al., 2019) and definitively qualifies the CFD to support the design of the future Gen IV reactors. One of the CRS4 tasks in the “Pool Thermal-hydraulics” work package of the SESAME project was to perform advanced CFD simulations of the liquid lead circulation under freezing condition, in collaboration with the Research Centre Rez (CVR). Therefore, for the solidification model development task, a dedicated numerical model representing the experimental facility used for the solidification experiments campaigns (Dofek et al., 2016) has been built. The computational model is then validated on the experimental data. The full 3D numerical model of the SESAME-Stand pool has been built in STARCCM+ (STAR-CCM+, 2019), versions 13 and 14. The model includes a pool–type vessel equipped with four heating rods, with an internal obstacle supporting the natural convection and with multiple thermocouple probes. The external surface of the experimental vessel is cooled down by air with controllable mass–flow rate which allows steady–state and transient studies of the lead solidification under various conditions. The model was developed first as a pre-test model and the preliminary results obtained were described in (Melichar et al., 2018) and in (Iannone et al., 2017). In the post-test phase, the model is further analyzed and improved. The detailed description of the model as well as the results obtained are included in (Profir and Moreau, 2019). Steady state and fast transient cases are reproduced according to the set of experimental tests performed at CVR. The thermal behavior for different heat sources and different air velocities was simulated using
Temperature gradient in the Argon layer Length of Argon layer Argon thermal conductivity [W/(m-K)] Added argon thermal conductivity [W/(m-K)] Heat transfer coefficient [W/m2-K] Heat flux [W/m2] specific heat, heat capacity [J/(kg-K)] radius [m] dynamic viscosity [Pa-s] thermal conductivity [W/(m-K)] density [kg/m3] turbulent viscosity [Pa-s] kinematic viscosity [m2-s]
input from the experimental work (Melichar et al., 2018). The pre-test simulations were partially satisfactory. An investigation campaign was actuated in order to understand the numerical-experimental differences and improve the numerical model. Initially, the numerical results presented large differences with respect to the experimental ones. The first experimental results obtained in the SESAME-Stand facility showed that the heat losses were much more intense than expected from the simulation and source of these additional losses should be detected. Observing that the air channel external border was much hotter than bare hand could stand, we understood that the important differences between the numerical and the experimental results could be attributed to thermal radiation. An important contribution was given by the radiation thermal flux, not initially modelled. Progressively improving the thermal radiation modelling and optimizing the mesh refinement, the differences reduced significantly and a satisfying agreement with the experimental resultswas obtained. The paper is organized as follows. After the initial preparation (geometry, mesh and physics) described in Section 2, the first pre-test results are briefly illustrated in Section 3. The first simulations performed with the experimental inputs and the consequent analysis of the results and of the differences against the experiment are described in Section 4, with particular attention to the thermal radiation issue. A set of sensitivity studies investigating the dependence on time step, mesh, view factor and Prandtl number are described in Section 5. The final numerical results and the comparisons with the experimental tests, illustrating the main features of each simulated case, in relation to the thermal radiation modeling, are illustrated in Section 6. Finally, the conclusions of this numerical and experimental campaign are summarized in Section 7. 2. SESAME-Stand facility and numerical model 2.1. Experimental test section The main part of the facility is the experimental vessel (EV) (Melichar et al., 2019). The internal diameter of the EV is 300 mm and the total height is 404 mm. The lead level is approx. 30 mm from the lid and the space over the lead volume is filled with inert gas (argon) to compensate volumetric changes of the lead. The bottom of the EV has thickness of 25 mm and is equipped with 2° chamfer to facilitate draining. The external diameter of the EV is 360 mm and the external surface is formed by 48 triangular ribs of height 19.6 mm and inner angle 60° for enhancement of the heat transfer area. The EV is made of cast stainless steel 316. A replaceable cylindrical obstacle made of stainless steel has height of 180 mm and is placed in the center of the EV. The internal diameter of the obstacle is 66 mm and the external one is 97 mm. Four main electric heating rods – main heaters (MH) are located inside the obstacle. Thermal power of each heater is controllable between 0 and 2 kW (total heating power is 8 kW). The total 2
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
length of a main heater is 335 mm, with heated length of 280 mm. 19 thermocouple wells with external diameter of 5 mm are distributed all over the EV (Fig. 1, left). The air channel (AC) surrounds the EV and is not insulated. There is no thermal instrumentation on the air channel walls. The air flow is driven by a radial fan placed below the vessel. The maximum air mass–flow rate is 0.5 kg/s. In the vertical position corresponding to the EV, the air channel is cylindrical and its external diameter is 450 mm. A conical cover protecting cables coming out of the EV is located in the AC under the EV. Another purpose of the cone is to ensure better air flow distribution. The upper part of the AC acts as a chimney providing ventilation of the hot air from the stand to the experimental hall. This part is equipped with a mass–flow meter (Fig. 1, right).
• • •
2.2. Computational model of SESAME-Stand pool in STAR-CCM+
2.2.2. Mesh Once the geometrical parts have been assigned to the regions and the interfaces between the various solids-fluids, solid–solid, fluid–fluid regions have been defined, the surface and volume mesh have been generated. In STAR-CCM+, the mesh operations have been transferred from regions to the geometrical parts levels. Here, different mesh operations can be defined for different geometry parts. The interfaces assure the continuity between the adjacent regions, although not conserving the mesh conformity when different regions are meshed with different operations. The mesh models used in the simulation are the Surface Remesher and the Polyhedral/Trimmed Mesher. The base size of the mesh was initially 4 mm and successively it was halved in the lead domain and the TC/heater pipes. The Polyhedral mesh model was used mostly in the pre-test simulations but, since it produced systematically crushes of the simulations at the beginning of the solidification, the Trimmed mesh model was adopted definitively. For the fluid regions, the Prism Layer Mesher is added with three prism layers with total thickness of 2 mm. The maximum Wall Y + values are in the range of 5–13. Taking advantage of the central symmetry, a quarter of the domain with suited symmetry plane boundary conditions, was considered in order to reduce the computational run time. This is to focus the
2.2.1. CAD geometry The 3D geometrical model was built with the STAR-CCM + 3D-CAD modeler, according to the real dimensions indicated in the experimental facility description (Dofek et al., 2016). The CFD model, illustrated in Fig. 2, includes the following components:
• Experimental vessel in which the thermal-hydraulics and freezing
• • •
distributed in the pool and one central TCs probe. They have a 6 mm diameter and a length of 200 mm. In the numerical model, the pipes hosting the thermocouples have a hexagonal geometry instead of circular, for better mesh matching at the interfaces with the fluid. Obstacle. An internal cylindrical obstacle is inserted in the vessel to enhance the natural circulation of the liquid lead. Lead. The lead domain is a cylinder with the diameter of 300 mm and the height 300 mm. The fluid volume is obtained by subtracting all the solids from the cylinder. The vessel’s internal is occupied by lead and argon. Argon layer. The argon volume occupies the upper part of the vessel’s internal for 30 mm of height above the lead volume. In the experiment, the argon layer is meant to compensate the volumetric expansion of the lead and to avoid air penetration in the vessel.
process take place. The vessel solid walls are fins in the form of triangular ribs for larger heat transfer area. The bottom portion of the EV has a thickness of 25 mm. The vessel main dimensions are: Inner diameter = 300 mm, Outer diameter = 360 mm, Height = 404 mm. Cooling air channel. The air channel main dimensions are: Inner diameter = 320 mm, Outer diameter = 450 mm, Height = 300 mm. Main heater formed by four heating rods placed in the center of the experimental vessel, around the central TCs probe. In the numerical model, the heating rods are simple steel cylinders with outer diameter 10 mm. The heater consists in a unique solid region with a volumetric heat source applied. Thermocouples probes. There are 12 bottom and 8 top TCs probes
Fig. 1. Layout of the thermocouples (left) and air channel (right). 3
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 2. CFD geometry. Top left: full CFD model. Top right: solid vessel. Bottom left: air channel, inlet in green, outlet in yellow. Bottom center: central and top TC pipes in violet, bottom TC pipes in orange, obstacle in dark violet, heaters in green, vessel in gray. Bottom right: fluid domains, Lead/Argon interface in pale gray. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 3. Volume mesh on a transversal section in the polyhedral model (right) and in trimmed model (right).
4
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
enthalpy. The steps of the algorithm were: i. Impose the cooling air velocity and temperature ii. First guess of the heat source based on the heat flux removed by the air iii. Reach the steady state solution iv. Look at the minimum temperature of the lead v. Adapt dynamically the heat source to bring the global enthalpy content to the value that provides the correct minimum temperature:
Table 1 Number of volume cells in each domain.
Air Lead Argon Solids Total
Polyhedral mesh
Trimmed mesh
Trimmed refined mesh
0.26 0.27 0.03 0.14 0.70
0.24 0.15 0.02 0.26 0.67
0.24 0.84 0.02 0.37 1.47
M M M M M
M M M M M
M M M M M
HS = min(20 kW, attention on the physics. In all the preliminary simulations with the whole domain we couldn’t notice asymmetric significant deviations. The total volume mesh representation consists of about 0.7 million of control volumes and about 1.5 million cells in the trimmed refined mesh case. Illustrations of the volume mesh are given in Fig. 3 while the number of cells for each domain is provided in Table 1.
max(HF
0. 02 × (H
H 0)), 0. 5 kW),
where HS is the heat source, HF is the heat flux removed by the cooling air, H is the total enthalpy of the lead, H0 is the target enthalpy of the lead. vi. Fix the new found heat source to a constant value. Once the global thermal equilibrium was reached, the thermal behavior for different heat sources and different air velocities was simulated. The numerical pre-test simulations showed a strange spotted behaviour in the very small velocity regions but, since the main flow was not perturbed, this issue was initially neglected. In fact, for small values, in the order of 3 cm/s, the velocity field couldn’t find a reasonable configuration (Fig. 4, left). Moreover, several technical issues were faced: systematic divergence leading to occasional and unexpected simulation crashes occurred when the solidification process started. Making a systematic analysis on a simplified (still 3D) case, some causes could be excluded: turbulence model, convergence parameters, mesh size, CFL and a very undue sensitivity to the reference density. Based on the Fluent model (Melichar et al., 2018); which has used a structured mesh to avoid simulation crash, the Polyhedral Mesh model was replaced with the Trimmed Mesh model which is the closer STARCCM + model to the structured mesh (Fig. 4, center). Instantaneously, the flow became defined (Fig. 4, right) and the solidification was well modelled with the default Melting-Solidification (sub-model of the Volume of Fluid framework) main settings. In STAR-CCM+, the access to the Melting-Solidification model is enabled once the Volume of Fluid (VOF) model is enabled in an Eulerian multiphase configuration with Lead as single phase. The cover gas is treated as a separated fluid continuum connected with the lead continuum by means of an impermeable baffle interface. In the numerical model, the volume fraction of lead occupies the full closed domain where the natural circulation and the freezing are simulated. The time model compatible with the VOF framework is the Unsteady State model. Only transient simulations are therefore performed. For the freezing simulation, the following numerical settings of the Melting-Solidification model are specified:
2.2.3. Physics The physical properties and the suitable numerical models are chosen. The following physical models are set up in the numerical simulation for the three fluid domains (lead, argon and air) and for the solid domains:
• Time: Unsteady • Material: Gas/Liquid • Flow: Segregated Flow • Equation of state: temperature dependent density for lead and argon, constant density for air • Viscous regime: Turbulent; Reynolds-Averaged Navier-Stokes • Reynolds-Averaged Turbulence: Realizable K-Epsilon model, TwoLayer All y + Wall Treatment • Eulerian Multiphase for the lead domain with Volume of Fluid (VOF) and Melting-Solidification models • Segregated Fluid Temperature model (as the Melting-Solidification model requires a non-isothermal energy model) • Constant density and Segregated Solid Energy models for the solid domains.
The physical properties of the fluids and solids involved are set up according to the recommendations in (Niceno and Badillo, 2017) and (Handbook, 2015). The working fluid in the vessel is the molten lead. The physical properties of the liquid lead are listed in Table 2. The thermo-physical properties of the lead in solid state are given in Table 3. The physical properties of the Argon are listed in Table 4. The vessel walls, the thermocouples pipes and the heater are made of stainless steel with the physical properties listed in Table 5.
• Latent heat of fusion = 23070/5 J/kg (20%) • Liquidus temperature = 600.6 K • Solidus temperature = 600.6 K • Flow stop relative solid fraction = 0.02 • Turbulent Prandtl number set to 10 where the flow is stopped
3. Pre-test numerical simulations Before proceeding with the effective freezing, it is convenient to bring the system close to the freezing onset in a state of thermal equilibrium. However, the steady state minimum lead temperature is very sensitive to the heat source value and the air velocity. We couldn’t find this value with trial and error tests and thus we used a dynamic algorithm based on the heat flux removed by the cooling air and on a target
6
(Flow Stop Flag greater than 0.5) and the default turbulent Prandtl number = 0.9 everywhere else.
In practice, the default settings of the VOF environment are kept.
Table 2 Liquid lead properties. Property
SI unit
Correlation
Melting pointTM Density Heat capacity at constant pressure
K kg/m3 J/kg-K
TM = 600.6 = 11441
1.2792T
cP = 176.2
4.923 × 10 2T + 1.544 × 10 5T2
Dynamic viscosity
Pa-s
µ = 4.55 × 10 4exp(1069/T) k = 9.2 + 0.011T
Thermal conductivity
W/m-K
5
Temp. range (K)
1.524 × 106T
2
TM TM
1900 1300
TM
1500
TM
1400
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Table 3 Solid lead properties. Property
SI unit
Heat capacity at constant pressure
J/kg-K
Thermal conductivity
W/m-K
Correlation
Temp. range (K)
cP = 126.49 + 4 .7× 10 k = 30
17.3 × 10 3 (T
2 (T
300
273.15)
600.6) = 40.39
17.3 × 10 3T
80
TM TM
Table 4 Physical properties of Argon at 1 bar. Property
SI unit
Correlation
Density
Kg/m3
Heat capacity Thermal conductivity
J/kg-K W/m-K
= 4.5751 nothing cP = 521
Dynamic viscosity
Pa-s
Temp. range (K)
1.57317 × 10 2T + 2.32033 × 10-5T2
1.24240 × 10-8T3
k = 1.55271 × 10-3 + 5.96287 × 10-5T -1.92121 ×10-8T2 µ=1.82652x10-6 + 7.69524x10-8T- 2.51082x10-11T2
SI unit
Correlation
Density Heat capacity at constant pressure Thermal conductivity
Kg/m3 J/kg-K
= 8090 cP = 500
0
W/m-K
k = 14.5 + 0.015T
100
700
273 273
700 700
273
700
difference compared with the experimental results, the temperature values being more than 200° higher than the experimental values (Fig. 6, right). Similar numerical results and similar discrepancies against the experiment were obtained also in the CVR model (Fig. 6).
Table 5 Physical properties of stainless steel. Property
273
Temp. range (C)
4. Analysis of the numerical-experimental differences
100
The experimental results showed that the facility heat losses were much more intense than expected from the simulation and we should look for potential origin of these added losses. Then, observing that the air channel external border was much hotter than bare hand could stand, it was observed that the important differences between the numerical and the experimental results could be caused by the thermal radiation influence, which was not yet contemplated in the model. It turns out that thermal radiation definitively plays an important role in the facility cooling and should absolutely be taken in consideration in its modelling. However, the Thermal Radiation model included in STARCCM + was not employed for some reasons:
500
The Liquidus and the Solidus temperature are both 600.6 K. Since the objective is to reach a steady state frozen fraction, for the acceleration of the solidification, the latent heat of fusion is set at 20% of the physical value. It must be stressed that, the liquid being a pure substance, the transient region is sharp and there is no mushy zone (Moreau et al., 2019). Moreover, there is no mesh adaptation used. In particular, no refinement near the front is dynamically performed. However, it could be a valuable direction for future improvement. In view of the comparison with the results of the experimental tests, the simulation is run with the following experimental input: Heat source HS = 5600 W and Air Mass Flow Rate MFR = 0.145 kg/s. The temperature field on the symmetry planes (in the lead and the vessel) and the maximum temperature curve approaching the steady state condition are shown in Fig. 5. The numerical results present a big
• the incompatibility with the geometry of our model, which consists •
in different types of surfaces and different domains (lead and cover gas); the uncertainty about the material emissivity, which has a large variation, depending on the treatment of the material (Table of
Fig. 4. Velocity field in the liquid lead domain in pre-test simulations with polyhedral mesh showing spotted behavior (left) and in simulations with trimmed mesh showing “clean” configuration (center and right). 6
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 5. Temperature field in the STAR-CCM + model with about 200° higher than experimental values.
Fig. 6. Temperature field in the Fluent model (left) and in the experiment (right).
Emissivity, 2003);
cover gas and at the external lateral walls of the vessel in contact with the air channel (Fig. 7). The thermal radiation is modeled as radiative heat fluxes at the vessel’s walls and as an increased thermal conductivity in the cover gas. More precisely, it is defined as follows:
• the thermal radiation that we modelled in an approximated manner
is sufficient for our model which at this stage is an acceptable compromise between the level of geometrical details, physical fidelity and computational costs.
i. Vessel’s bottom radiative heat flux o φ = -e* σ (T4-T4bottom), where T is the local temperature in Kelvin and: ▪ e*=1/(1 + 2(1-e)/e) = 3/5 = 0.6 is the reduced emissivity corresponding to the material emissivity e = 0.75 ▪ σ = 5.67x10-8 W/m2-K4 is the Stefan-Boltzman constant ▪ Tbottom = aT+(1-a)Tair-In (a = 0.5), Tair-In = 20 °C ii. Vessel‘s top radiative heat flux o φ = -e* σ (T4-T4top), where
4.1. Thermal radiation modelling The thermal radiation issue was already addressed in (Moreau et al., 2019), in the CIRCE facility modeling, where the feature was implemented as an increased temperature dependent thermal conductivity. The techniques developed for CIRCE were adapted to the SESAME-Stand model. Thermal radiation is therefore defined and implemented at the top and at the bottom of the steel vessel, in the band of 7
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 7. Boundaries and regions where radiative heat fluxes are defined.
Fig. 8. Temperature field on symmetry planes in Case I (left) and Case II (right).
▪ Ttop = bT+(1-b)Tair-Out (b = 0.5), Tair-Out = 50 °C iii. Lateral radiative heat flux applied at the external wall of the vessel in contact with the air channel o φ= - α e* σ (T4-T4EW), where ▪ TEW = 70 °C is the temperature on the external surface of the air channel provided by the experimental measurements and =1/2 is the view factor of the opposite face, to take into account the triangular ribs shape of the external vessel’s surface iv. Heat flux approximated as an increased thermal conductivity in the cover gas o φ = e* σ (T14-T04)≃ 4 e* σ T3 (T1-T0)≃4 e* σ T3 l T n, where l = 0.03 m is the length of the argon layer, T is the temperature gradient in the argon layer and n is normal direction o the argon thermal conductivity k0 = 0.04 W/m-K is increased by k+= 4e*σT3 l W/m-K.
described in Section 4.1), the simulations with different inputs from the experimental tests could be run. The first pre-test simulation was continued with the thermal radiation modelled. A second case, more representative of the freezing front propagation, was also simulated with the thermal radiation included. All experimental measurements have shown that is reasonable to consider that the results are compliant with the geometrical symmetries. For this reason, the comparison between the experimental TCs data on a plane with the numerical TCs data on two half planes is meaningful. The numerical results obtained in these two cases, together with their main characteristics, are illustrated hereafter. Case I.HS = 5.6 kW and MFR = 0.145 kg/s
• Maximum temperature: 30 K above the experimental value (Fig. 8, Fig. 9) • Minimum temperature: 330 °C, in agreement with the experimental values • Numerical solidification fraction: 0.22%. Experimental frozen frac-
The model is then run with the thermal radiation included. The numerical results show better agreement with the experimental results. Including the radiations, the temperature profile becomes much more similar and the difference of temperature falls considerably.
tion: 5%
Case II. HS = 4 kW and MFR = 0.101 kg/s
• Maximum temperature: 5 K below the experimental value at the
4.2. Improved numerical simulations based on experimental input
center domain and 30 K above the experimental values near the wall
After having gradually improved the heat transfer modelling (as 8
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 9. Temperature on representative TCs in Case I. Numerical values (left) and experimental values (right).
Fig. 10. Temperature on representative TCs in Case II. Numerical values (left) and experimental values (right).
(Fig. 8, Fig. 10)
partial frozen fraction of 12.5% is well suitable for these purposes.
tion: 45%
5.1. Time step dependence study
• Minimum temperature: 30 K above the experimental values • Numerical solidification fraction: 12.5%. Experimental frozen frac5. Sensitivity/optimization studies on Case II. HS = 4 kW
First of all, a time step dependence study was performed on the steady state Case HS = 4 kW comparing the MFR curves and solidification percentage curves with different time steps. The study revealed a certain dependence of the time step; the solidification is slightly faster with halved time step and the curve of solidification is correlated with the variation of the mass flow rate, as can be noticed in the curves in Fig. 11.
Since the first steady state cases with the thermal radiation implemented gave results still different from the experiment, Case II. HS = 4 kW was investigated further. Following the engineering approach of this numerical work, a series of sensitivity and optimization studies were performed, looking at the effect of parameters as time step, mesh size, material emissivity or view factor, this last one having a noticeable influence in the thermal radiation modelling. The case with 9
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 11. Lead MFR (left) and solidification fraction (right) with TS = 1 s and TS = 0.5 s in the Case HS = 4 kW.
Fig. 12. Lead MFR (left) and solidification fractions (right) with TS = 0.5 s and initial and refined mesh.
Fig. 13. Solidification fractions comparison with different mesh and different turbulent Prandtl number: Prt = 0.9 until time 360 min followed by Prt = 2, Prt = 10 and Prandlt number in Kays correlation.
5.2. Mesh dependence study
the volume mesh increased from 0.7 M to 3.5 M cells and the solidification fraction increased from 12% to 19%. The difference is therefore noticeable, about 56%, as shown in Fig. 12 (right). This is an indication that for the future cases we need to find a compromise between the
The next sensitivity study concerned the mesh dependence solution. Dividing the base size of the surface mesh by two (from 4 mm to 2 mm), 10
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 14. Definition of the view factor on the vessels ribs. Table 6 Comparative results with different view factor and different mesh. Case HS = 4 kW
Max Lead Temp (C) Frozen fraction % Lead MFR (kg/s)
View factor with Initial Mesh
View factor with Fine Mesh
α = 1/2
α = f(r)
α = 1/2
α = f(r)
364
361
363
360
371
12
15
19
22
45
0.186
0.192
0.188
0.2
Table 8 Input and results in the numerical steady state cases compared with experimental results.
Experimental results
Case
1 2 3 4 5
Table 7 Numerical steady state and transient cases.
Steady State Cases 1 2 3 4 5 Transient Cases 1 2
Heat source (kW)
Air MFR (kg/s)
5.6 5.6 4.96 4 5.6
0.145 0.167 0.101 0.101 0.095
0.4 5.6
0.095 0.095
Without Top Lateral heat transfer
With Top Lateral heat transfer
Experiment
HS (kW)
Air MFR (kg/s)
HTC (W/ m2-K)
Frozen fraction %
HTC (W/ m2-K)
Frozen fraction %
Estimated Frozen fraction
5.6 5.6 4.96 4 5.6
0.145 0.167 0.101 0.101 0.095
0 0 0 0 0
1.5 6 0 34 0
33 36 6.7 13 12
4.4 12 0.2 58 0
5% 16% 0 45% 0
investigated. It turns out that un-refining the air channel does not bring any change to the lead temperature. This was an indication that we could refine the mesh only in the lead domain and reduce the computational cost considerably. Therefore, the simulation proceeded with the refinement of the lead domain only and the mesh independence of the freezing curve was confirmed. Indeed, refining the mesh when the freezing was at 12% and continuing the simulation, the previously obtained frozen fraction of 19% with globally refined mesh (with 3.5 M cells) was obtained with 1.4 M cells, with an important saving of computational cost. It is recognized that it’s not known whether reasonable mesh independence is already reached. However, we are aware that mesh independence in the lead domain could not be demonstrated.
computational cost and a reliable convergence. The Lead MFR is sensitive to the mesh refinement as well, being slightly higher in the initial stage of the freezing and almost the same at the steady state, as indicated in Fig. 12 (left). The higher frozen fraction percentage with the finer mesh is consistent with the lower values of the minimum temperature obtained with the finer mesh. The mesh dependence is a serious issue because the finer mesh requirement leads to 3.5 M control volumes for only a quarter of the domain simulated. This means that an eventual upgrade to the entire domain to study the future asymmetric case becomes challenging. As most of the control volumes are in the air channel, separate effect of the lead domain (and partially of the solid domains) mesh density was
5.3. Sensitivity to Prandtl number Finally, a sensitivity study was performed by varying the value of the turbulent Prandtl number Prt. In the solidification model, the turbulent quantities are “frozen” together with the flow stop criterion and the solid lead keeps some turbulent part of the conductivity. This is numerically corrected by dynamically setting the turbulent Prandtl number at 106 where the flow is stopped. The effect is to reduce the turbulent convection part of the heat transfer to a negligible amount. In 11
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 15. Temperature field on symmetry planes in Case 1 (left) and Case 2 (right).
Fig. 16. Case 1. Final temperatures on representative TCs numerical (left) and experimental (right).
the sensitivity study performed on the Case HS = 4 kW, starting from the time 360 min, when the steady state was essentially reached, the turbulent Prandtl number was increased from the initial 0.9 value first to 2 and then to 10. The solidification fraction increases noticeably with the increment of the turbulent Prandtl number, both in the case with the initial mesh and with the refined mesh, as shown in Fig. 13. We observe that the use of different Prandtl number values induces a large influence on the results. In order to decide which the more correct value is, we need to compare with more elaborated Prandtl number models. The most elaborated and precise formula for the turbulent Prandtl number making use only of local variables, and thus suited for CFD implementation, is the Kays correlation (Bartosiewicz et al., 2014); (Kays, 1994): Prt = 2/Pet + 0.85, where Pet is called the turbulent Peclet number and is function of the Prandtl number and of the ratio of turbulent viscosity t and kinematic viscosity : Pet = Pr t / = t cp/ k . Being built on theoretical considerations, the Kays correlation is expected to be more reliable than a mere constant value. Therefore, in
absence of experimental data, we check the validity of different constant values implementation with the Kays correlation. The results obtained with the Kays correlation fully confirm the results obtained with Prt = 0.9, as illustrated by the solidification curves in Fig. 13, which show how almost nothing changes from the moment 360 min when the turbulent Prandtl number is modified. 5.4. Sensitivity to thermal radiation modelling Initially, in the definition of the radiative heat fluxes, the view factor of the vessel’s wall surface was set to α = 1/2 to take into account the orientation of the normal. After that, the view factor was defined as a function of the depth of the vessel’s wall in the interval [1/ 3,2/3], as illustrated in Fig. 14. The results obtained with this definition of the view factor are given in Table 6. A small increment of the solidification fraction was registered, both with the initial mesh and with the refined mesh, getting closer to the experimental results. 5.4.1. Top lateral heat transfer improvement During this investigation, it was revealed that the effect of the view 12
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 17. Case 2. Final temperatures on representative TCs numerical (left) and experimental (right).
Fig. 18. Evolution of the frozen fraction in Case 4 from 0% to 58%: 12% with initial mesh, 20% with refined mesh, 34% with adjusted view factor, 58% with improved thermal radiation.
factor is quite important and the focus was put on the correct implementation of this factor. Initially, as described in Section 4.1, the view factor used in the definition of the radiative heat flux was α = 1/ 2, both on the triangular ribs side wall and on the slanted portion of wall at the top in contact with the air. Erroneously, in the slanted portion, the view factor was α = 1/2; moreover, at the top vertical surface not in contact with the air channel (the vertical top surface in Fig. 7), the thermal radiation was completely omitted. Therefore, the thermal radiation modelling had to be further improved. The view factor in the lateral ribs is set to vary linearly with the radius, α = f(r) in the interval [1/3,2/3]. It is kept constant at α = 1 in the slanted and vertical top portions of the vessel’s wall. In the second step, the heat flux at the top lateral wall of the vessel, above the air channel (the top vertical portion indicated in Fig. 7), was improved. Additionally to the already applied radiative heat flux defined according to the formula φ = -αe*σ(T4-TEW4), α = 1, a heat flux component based on the heat transfer coefficient of the slanted portion,
calculated in the following manner, was added: i. the heat flux q through the top slanting portion of wall in contact with the air channel was evaluated from the measured heat transfer through the surface and the surface area; ii. the heat transfer coefficient of this portion was obtained, h = q/ΔT, by dividing the heat flux to the bulk temperature ΔT = Tmean-Tairout, where Tmean is the average temperature on the slanting wall and Tair-Out = 50 °C is the air outlet temperature; iii. the obtained heat transfer coefficient is used in the calculation of the heat flux added to the top wall: φ = -h(T-Tair-Out), where T is the local wall temperature in Kelvin. The method to define the convective heat transfer coefficient (HTC) is very crude, since it implicitly relies on a fully developed flow assumption, which is not our case. It is to be considered as an initial guess. The results obtained from the five cases and the sensitivity study 13
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 19. Temperature field on symmetry planes in Case 3 (left) and Case 4 (right).
Fig. 20. Case 3. Final temperatures on representative TCs numerical (left) and experimental (right).
on the HTC confirm that the adopted method is at least consistent. The HTC was calculated from the heat transfer including both the radiative and conductive components of the slanted wall, while only the conductive component should have been considered. However, the change of shape of the channel is very likely to induce a strong local turbulence increase and consequential increase of the HTC. Therefore, there is probably a compensation of errors beneficial to the determination of the heat transfer coefficient.
partially refined mesh and adjusted view factors and ii) top thermal radiation modelling, the five steady state cases and the two transient cases in Table 7 were simulated. 6.1. Numerical steady state results In this section, the main results of the five steady state cases simulated are illustrated and the comparison with the experimental results is realized. The percentage of frozen fraction obtained in each case is anticipated in Table 8 in comparison with the frozen fraction obtained in the experiment. The frozen fractions from the experiment are only indicative, being estimated from the measured TC values (Melichar et al., 2019). The main features of the simulated cases are illustrated. The temperature field on symmetry planes and the temperature distribution are
6. Numerical results. Comparison with the experimental results Based on the gradually improved preliminary post-test simulation as described in Section 4.1 and Section 5.4.1, additional simulations following a part of the experimental matrix, with steady state cases and transient cases, were run (Melichar et al., 2019). Therefore, with i)
14
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 21. Case 4. Final temperatures on representative TCs numerical (left) and experimental (right).
Fig. 22. Case 5. Final steady state temperatures (left) against the experiment (right).
• Maximum temperature: in agreement with the experimental values (Fig. 15, Fig. 17) • Solidification fraction: 12% with implemented thermal radiation,
shown for each case in the pictures from Fig. 15 to Fig. 23. The thermal radiation model used to produce these pictures is the improved radiation model defined in Section 5.4.1. Case 1. HS = 5.6 kW and MFR = 0.145 kg/s
16% estimated experimentally
• Maximum temperature:
Case 3. HS = 4.96 kW and MFR = 0.101 kg/s
o 10 K above the experimental values in the case with still incomplete thermal radiation o in agreement with experimental values in case with improved thermal radiation (Fig. 15, Fig. 16) Solidification fraction: 4.4% with implemented thermal radiation, 5% estimated experimentally
• Maximum temperature:
o 15 K above the experimental values in the case with still incomplete thermal radiation o in agreement with the experimental value with improved thermal radiation (Fig. 19, Fig. 20) Solidification fraction: 0.2% numerically, 0% experimentally
•
•
Case 2. HS = 5.6 kW and MFR = 0.167 kg/s
15
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 23. Case 5. Final steady state temperature and velocity.
Case 4. HS = 4 kW and MFR = 0.101 kg/s
Nevertheless, the temperature curves cannot give a deep insight on the evolution and behavior of the freezing front. One can get this insight by looking at the pictures captured in some key moments which clearly give indications of the asymmetry of the solidification and melting phenomena, as illustrated in Fig. 26, Fig. 27 and Fig. 28.
• Maximum temperature: in agreement with the experimental value (Fig. 19, Fig. 21) • Solidification fraction (Fig. 18):
o 12% with initial mesh and 20% with refined mesh, without thermal radiation o 34% with adjusted view factor (but still incomplete thermal radiation) o 58% with improved thermal radiation, 45% estimated experimentally
7. Discussions and conclusions A computational model of the SESAME–Stand experimental facility pool, able to simulate the behavior in case of solidification of the liquid lead, was developed using the STAR-CCM + simulation code. The 3D simulations in buoyancy driven flow proved to be initially problematic with occasional and non-averted simulation crashes when the solidification process started. Proceeding with a systematic analysis, many causes could be excluded. Taking inspiration from the CVR model experience, the Trimmed Mesh model was used, and the technical issues have been overcome. The model became thus suitable for the simulation of the lead solidification experimental tests performed at CVR and the numerical implementation of the experimental matrix could be performed. The comparison between the first post-test numerical results and experimental tests results revealed initially strong differences (in the range of 200–300 K). Similar differences were also present in the CVR numerical model. Both partners agreed that the differences were probably to be mainly attributed to the fact that the thermal radiation was not taken into account. Thermal radiation was then modeled and included in the model as an increased temperature dependent thermal conductivity in the cover gas and by means of heat fluxes appropriately defined at the vessel’s walls (external top, bottom and lateral walls, cooling air channel interface walls). The heat losses at the bottom, top and side walls were estimated and variable view factors of the ribs surface have been implemented. The correct definition of the radiation came step by step, based on the experimental measurements on the air channel external wall. Once the thermal radiation was properly defined over the whole surface of the vessel, the differences between the numerical and experimental results became reasonable. Five steady state cases as well as one fast transient simulating the freezing followed by melting of the lead content based on the experimental matrix were simulated and compared with the experimental tests. The numerical steady states results show good agreement with the experimental results. The main features of each case are highlighted
Case 5. HS = 5.6 kW, initial condition for transient This case is the steady state used for the fast transient. It was obtained from the Case 1. HS = 5.6 kW, keeping the heat source constant and reducing the air mass flow rate from 0.145 kg/s to 0.095 kg/s. The case was first run with the initial mesh and with constant view factor. The comparison with the corresponding experimental test results shows quite similar trend of the temperature on the monitored thermocouples but with the following differences: Maximum temperature above the experiment by 25 K and Minimum temperature above the experiment by 10 K. With the refinement of the mesh and the correction on the view factor, the maximum temperature decreased by 15 K. With the further implementation of the heat flux at the top lateral wall, the temperatures are well aligned with the experimental values, as shown in Fig. 22. The final steady state temperature and velocity fields are illustrated in Fig. 23. The heat losses were measured, distinguished by surfaces and types of heat transfer; this distinction provided an indication on the percentage of the thermal radiation, which proved to be roughly 50%. 6.2. Numerical transient results The first transient case was simulated by suddenly decreasing the heating power in the previous case, from 5600 W to 400 W. The temperature of the lead decreases and almost the entire content of lead is solidified. The heating power was switched again to 5.6 kW and the lead temperature rose sharply as shown in Fig. 24. The comparison with the experimental fast transients shows quite similar behavior, as illustrated in Fig. 25. As the comparison reveals, the numerical results are in excellent agreement with the experimental temperatures.
16
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 24. Numerical steady state (1 h), followed by fast transient with solidification (< 1h) and re-melting (4 h).
Fig. 25. Experimental steady state (1 h), followed by fast transient with solidification (< 1h) and re-melting (4 h).
17
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al.
Fig. 26. Initial condition at time 1:04 (left), 55% freezing at time 1:30 (center), 73% freezing at time 1:38 (right).
Fig. 27. Melting kick-off 66% freezing at time 1:41 (left), 64% at time 1:43 (center), 44% at time 2:18 (right).
Fig. 28. Melting in final phase 26% at time 2:54 (left), 8% at time 3:54 (center), 4% at time 4:18 (right).
and the comparisons with the experimental results are illustrated by means of temperature plots and data for representative thermocouples. The obtained results confirmed that the thermal radiation has an important contribution since it is responsible for about 50% of the heat losses. Finally, the fast-transient experimental tests with sudden decreasing and increasing of the heat source were reproduced numerically. It can be concluded that the numerical simulations provided good agreement with the experimental results, especially if it is considered the high sensitivity of the frozen fraction to the heating source, convective heat transfer change with the cooling air velocity and change of thermal radiation heat transfer based on the change of the vessel outer walls temperatures.
Acknowledgments The work described in this paper has been funded by Autonomous Region of Sardinia and the EURATOM research and training program 2014-2018 under the grant No. 654935 (SESAME). References Bartosiewicz, Y., Duponcheel, M., Manconi, M., Winckelmans, G., Bricteux, L., 2014. Turbulence modelling at low Prandtl number, Chapter in the book: Fluid Mechanics and Chemistry for Safety Issues in HLM Nuclear Reactors. von Karman Institute for Fluid Dynamics, Belgium ISBN-13 978-2-87516-064-5. STAR-CCM+ User’s Guide Releases 13 and 14, Siemens. PLM Software. Dofek, I., Frybort, O., Vit, J., Slinc, M., 2016. “Description of the Meliloo – stand test facility and test matrix”, SESAME project deliverable D3.8, https://cordis.europa.eu/project/rcn/ 198041/results/en. Handbook on Lead-bismuth Eutectic Alloy and Lead Properties, 2015. Materials Compatibility, Thermal-hydraulics and Technologies, https://www.oecd-nea.org/science/pubs/2015/ 7268-lead-bismuth-2015.pdf. M. Iannone, I. Dofek, T. Melichar, W. Borreani, G. Lomonaco, V. Moreau, M. Profir, 2017. Development of CFD models and pre-test calculations for thermal-hydraulics and freezing experiments on Lead coolant, Proceedings of the International Conference Nuclear Energy for New Europe, Bled, Slovenia, September 11-14, 2017, http://publications.crs4.it/ pubdocs/2017/IDMBLMP17/.
CRediT authorship contribution statement Manuela Profir: Writing - original draft, Writing - review & editing, Investigation, Methodology, Visualization. Vincent Moreau: Conceptualization, Supervision, Project administration. Tomáš Melichar: Data curation, Validation. 18
Nuclear Engineering and Design 359 (2020) 110482
M. Profir, et al. Kays, W., 1994. Turbulent Prandtl number-where we are. J. Heat Transfer-Trans. ASME 116 (2), 284–295. Koloszar, L., Moreau, V., 2016. (U)RANS pool thermal hydraulics, Chapter in the book Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors, Woodhead Publishing Series in Energy. ISBN: 978-0-08-101980-1(print), 978-0-08-101981-8 (online). Melichar, T., Belka, M., Frybort, O., Kordac, M., Moreau, V., Profir, M., 2018. Meliloo: pre-test simulations and measurements, SESAME project deliverable D3.10, http://publications. crs4.it/pubdocs/2018/MBFKPM18a/. Melichar, T., Frybort, O., Matous, P., Belka, M., Moreau, V., Profir, M., 2019. Design and first experimental data from SESAME-Stand for lead solidification experiments, SESAME International Workshop, Petten, The Netherlands, 19-21 March, http://publications.crs4. it/pubdocs/2019/MFMBMP19/. Melichar, T., Matous, P., Achuthan, N., Belka, M., Frybort, O., 2019. FLUENT solidification simulations, SESAME project deliverable D3.15, https://cordis.europa.eu/project/rcn/ 198041/results/en. Moreau, V., Profir, M., Alemberti, A., Frignani, M., Merli, F., Belka, M., Frybort, O., Melichar, T., Tarantino, M., Franke, S., Class, A., Yanez, J., Grishchenko, D., Jeltsov, M., Kudinov, P., Roelofs, F., Zwijsen, K., Visser, D., Badillo, A., Niceno, B., Martelli, D., 2019. Pool CFD Modelling: Lessons from the SESAME Project, NED, Vol 355, 2019 doi: 10.1016/j.nucengdes.2019.110343. Moreau, V., Towards a full 3Dsimulation of MYRRHA primary coolant loop, pp. 1-30, chapter in the book: Fluid mechanics and chemistry for safety issues in HLM nuclear reactors, VKI,
2014, ISBN-13 978-2-87516-064-5. B. Niceno, A. Badillo, 2017. “Summary for the recommended physical properties for the main fluid and liquid materials used for the MELILOO-stand and TALL-3D freezing experiments.” SESAME project deliverable. Profir, M., Moreau, V., 2019. Solidification model development, SESAME project deliverable D3.13, https://cordis.europa.eu/project/rcn/198041/results/en. Roelofs, F., Shams, A., Moreau, V., Di Piazza, I., Gerschenfeld, A., Planquart, P., Tarantino, M., 2016. Liquid metal thermal hydraulics, state of the art and beyond: the SESAME project. In: ENC 2016, Warsaw, Poland. Roelofs, F., Gerschenfeld, A., Tarantino, M., Van Tichelen, K., Pointer, W.D., 2016, Thermalhydraulic challenges in liquid-metal-cooled-reactors, Chapter in the book Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors, Woodhead Publishing Series in Energy. ISBN: 978-0-08-101980-1(print), 978-0-08-101981-8 (online). Table of Emissivity of various surfaces, https://www-eng.lbl.gov/, 2013. SESAME project http://sesame-h2020.eu/, 2019. Visser, D., Profir, M., Moreau, V., 2019. CFD model of ALFRED primary loop, SESAME project deliverable D3.7. Zwijsen, K., Dovizio, D., Moreau, V., Roelofs, F., 2019. CFD modelling of the CIRCE facility. NED 353. https://doi.org/10.1016/j.nucengdes.2019.110277. Zwijsen, K. et al., 2018. Numerical Simulations at Different Scales for the CIRCE Facility, ICAPP 2018.
19