Numerical study of the effect of geometry on the behaviour of internally heated melt pools for in-vessel melt retention

In-vessel melt retention is a vital safety design that is incorporated in some nuclear reactors to protect against beyond-design basis accidents. This entails the retention of the molten fuel, colloquially known as corium, to be retained within the reactor pressure vessel (RPV) for extended periods of time without any release of radioactivity to the surroundings. Multiple experiments have been done to better understand this phenomenon, such as the COPRA experiment, which is simulated in this study. The study aims to analyse the COPRA experiment’s forced convection cooling (via water) while also simulating cooling via natural convection. After comparing this to the original COPRA experiment’s data, the study showcases the effect of geometry variation in such experiments by testing different RPV geometry sizes. Our results can better inform the design of future severe accident experiments while taking into account the historical aspects of severe accident research, as the popularity of Small Modular Reactors (SMRs) increases. Our analysis of the variation in geometry provides insights for RPVs that are significantly smaller than the traditional pressurised water reactors’ RPV size. The difference in cooling effect and crust formation is presented for this variation of vessel size. The difference is not significant but is noteworthy for the design of future experiments and reactors.


Severe Accident Management strategy (SAM)
The meltdown of the reactor core containing the fuel elements followed by the relocation of the core is a beyond design-basis severe accident scenario in reactors, due to which large pool of reactive materials, colloquially known as corium, accumulates in the lower plenum of the reactor pressure vessel (RPV). The main objective of nuclear safety is to prevent the release of corium, and its by-products, deposited at the bottom vessel, from escaping to the environment.
The Defense-in-Depth (DiD) approach covers beyond design basis accidents (BDBA) which structures the approach that should be taken by plant designers and operators to mitigate the risk of accidents in a facility. The severe accident management (SAM) strategy has evolved over the years to explore various safety avenues, enhanced by probabilistic safety assessments (PSAs) (Sehgal et al., 2011).
The DiD approach provides 3 physical barriers to prevent the release of radioactivity into the environment, which involves reactor design the corium in the lower head of the reactor pressure vessel (RPV) by flooding the whole lower plenum. EVMR, meanwhile, aims to collect and eject the corium into a dedicated sacrificial core catcher system. This core catcher is aimed to arrest, spread, and cool the corium. This study focuses on the IVMR strategy.
The IVMR strategy was originally developed as a way to enhance the safety of older PWRs in order to upgrade them to the enhanced safety standards. The Loviisa VVER-440 reactor was retro-fitted with this design in order to enhance its safety, in light of the Chernobyl accidents.
The IVMR strategy is encapsulated in a three-fold tiered approach as well, each with its own unique set of challenges. Ma et al. (2016) elaborate on the historical approaches of IVMR and state the general concepts that are still followed by this SAM strategy. They enumerate the approaches as: 1. in-situ core quenching, 2. particulate bed cooling, and 3. melt pool cooling The last approach of melt pool cooling is supposed to be activated in case all the systems fail and it aims to cool the molten corium pool in the lower plenum of the RPV.
The final approach of melt pool cooling aims to remove the decay of the melt pool with the help of a coolant flow. Succinctly speaking, this strategy aims to ensure that the heat generated in the pool should not exceed the cooling capacity of the vessel. The critical heat flux is of paramount importance here and governs these phenomena. Given the uncertainty involved in the heat generation, and subsequently the heat flux in the vessel due to the transient nature of this phenomenon, this simplistic sounding task proves to be an exceedingly difficult parameter to estimate. Sehgal et al. (2011) also describes how the formation of an oxide layer at the top of a melt pool can be a concern for vessel integrity due to the reduction in heat dissipation. Numerous experimental studies have been designed to explicitly study this phenomenon. The experiments use a range of simulant materials to better understand this phenomenon and better understand the thermofluid dynamics of this SAM system.

Internally heated melt pool experiments
Experiments such as COPO, ACOPO, SIMECO, LIVE, BALI, and RAS-PLAV were some of the many experiments conducted (Wei et al., 2018). Initially, these experiments were focused on yielding heat transfer correlations for demonstration in safety cases. For instance, the COPO experiment used water as a simulant and was devised specifically for the Loviisa Nuclear Power Plant (NPP) (Wei et al., 2018). The BALI experiment is notable for using a 1:1 RPV for the experiment (Bonnet, 1999) while the RASPLAV experiment is notable for extensively analysing the behaviour of melt pool stratification due to oxidation of zirconium (Asmolov et al., 2001;Sehgal et al., 2011;Aksenova et al., 1994;Strizhov, 1994;Behbahani et al., 1999).
Some experiments using molten salts consisting a mixture of KNO 3 and NaNO 3 have been employed as a substitute for corium in experiments such as LIVE (Miassoedov et al., 2014) and COPRA . The LIVE suite of experiments was conducted to investigate the natural convection behaviour of melt pool in a three-dimensional hemisphere where either eutectic or non-eutectic mixtures of KNO 3 −NaNO 3 at various molar ratios were adopted as the simulant material (Miassoedov et al., 2014). Similarly, The COPRA tests were conducted to study the natural convection behaviour of melt pool in a two-dimensional section of a quarter RPV slice where water tests and molten salt tests were performed in a sequence (Zhang et al., 2016a;Luo et al., 2018;Zhang et al., 2018).

COPRA experiment
The COPRA test vessel was a 2D quarter circular slice intended to replicate the RPV lower plenum of the HPR-1000 reactor (Zhang et al., 2016e,d,c,b). A schematic of the test setup is shown in Fig. 1. The internal Rayleigh number (Ra ′ ) value in the COPRA could potentially reach 10 17 , which corresponds to the prototypical reactor state. Based on the COPRA study, it is thought that the uncertainties of the IVMR research outcome can be decreased. This study is also expected to contribute to the theoretical understanding of transient and steady-state core melt behaviour in the reactor lower head, which is crucial for international IVMR research.
The vessel's inner radius and thickness were 2.2 m and 20 cm, respectively. The thermally insulated vertical wall was 25 mm thick, while the curved vessel wall was 30 mm thick. Twenty electrical heating rods were used to imitate decay heat, and these heating rods were spread horizontally across the tank. The electrical heating rods were used to heat the pool during the experiment. Direct heat convection between the pool and the vessel wall section below the pool surface transferred heat from the pool to the vessel wall. Furthermore, the heat was dissipated into the environment via a free surface at the top. The coolant outside the vessel wall chilled the vessel wall even further, when employed .
Water was used as the initial simulant in the first set of experiments (Zhang et al., 2016b). The data collected included melt pool temperature distribution, heat flux distribution along the vessel wall, and heat transfer towards the curved wall. Thermal stratification was observed in the melt pool at the quasi-steady-state and the temperature in the upper section of the melt pool was homogeneous. The heat flux increased as the polar angles rose along the curved wall. Heating powers and pool heights had no significant effect on the normalised temperature (ratio of local temperature to maximum temperature) and heat flux distribution, indicating that such experimental relations might be better suited for higher power density and other vessel dimensions. This series of tests yielded temperature and heat flow distribution relations, as well as − ′ (Nusselt and modified Rayleigh number) correlations.
Eventually, a mixture of 20 mol% NaNO 3 and 80 mol% KNO 3 (noneutectic salt) as well as 50 mol% NaNO 3 and 50 mol% KNO 3 (eutectic salt) were used as simulant material in COPRA salt tests. Similar investigations from the water test as well as other large-scale melt pool natural convection properties were explored in the experiment. The central melt pool, cooling channel, and higher cover were all part of the test facility. Two types of lids were employed to replicate top cooling and insulation conditions. The COPRA salt test had a similar heat transfer performance to the COPRA water test. For the salt test matrix, a solid crust was observed along the vessel wall. The melt pool temperature in the COPRA molten salt test was significantly higher than what was observed in the CO-PRA water test. The difference in temperature between the maximum pool temperature and the melt/crust interface temperature was used to calculate the Nusselt number. The impacts of upper cooling, melt material relocation position, melt pool height, and internal heating power level were explored in the COPRA molten salt test. During the quasi-steady-state, thermal stratification was also observed. The effect of non-central and central pouring was studied as well in the test matrix. The experimental results were also compared to those of RASPLAV, ACOPO, and BALI. The test matrix also explored the change that heating power could have on the solidified salt crust.
There is plausible evidence that the geometry size and internal Rayleigh number have a significant effect on the natural convection, heat flux, and the velocity distribution, as stated by Zhang et al. (2015) and as demonstrated by Wei and Chen (2019b) for COPRA. Wei and Chen (2019b) conducted the geometric variation study for COPRA test using water as the simulant material. A total of 7 tests, with increasing geometric radii of the vessel, were performed by them.
A. Venkateshwaran et al. Three tests having a radius less than 1 m (0.1 m, 0.3 m, 0.5 m) were grouped as a set and the remaining four as another set (1.0 m, 1.5 m, 2.0 m, 2.5 m). The natural convection behaviour of the water pool was similar within a particular set of geometric sizes, however the effect was different when the two sets were compared. The velocity and temperature field of water pool and the natural convective behaviour of different internal Rayleigh numbers was also studied by Wei and Chen (2019b). For all the geometric sizes, the internal Rayleigh number was found to be between 10 11 (for the smallest vessel) and 10 18 (for the largest vessel). This paper follows a similar investigation of the geometric size of RPV. Wei and Chen (2019b) use water, which is the most widely and convenient simulant, however the formation of crust at a high temperature and the tendency to model the real-world scenario has the potential to be better if salt is used.
This study focuses on 50 mol% NaNO 3 and 50 mol% KNO 3 (eutectic salt) used in COPRA as the simulant material and provides additional data for different radii of vessel sizes. We observe a dependency that the internal Rayleigh number and the geometrical size of the vessel exert on the convective behaviour of the melt pool. We intend to further the findings of the COPRA study which deals with a similar hypothesis Chen, 2020, 2019b) and explore the characteristics of the same with eutectic salt as the simulant material.
In this paper, − SST model for turbulence characteristics (in ANSYS Fluent) along with the solidification/melting models for crust formation characteristics are employed. It is first compared with the COPRA eutectic salt test's experimental data by using the commercial software ANSYS FLUENT (Ansys ® , 2021). Subsequently, the validated model is extended to study the geometric effect of varying vessel sizes (0.5 m, 1 m, and 2.2 m) with the same eutectic salt mixture of 50 mol% KNO 3 -50 mol% NaNO 3 as the simulant.

Simulation details
The natural convection behaviour of an internally heated melt pool is a complex behaviour that involves heat convection mechanism, fluid circulation, crust development, and phase change behaviour to mention a few.
In this study, internally heated melt pool natural convection behaviour is modelled using ANSYS FLUENT 21.0 (Ansys ® , 2021). Unlike the Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES), which requires a significant amount of processing resources, RANS turbulence models are computationally cheaper.
Wei and Chen (2019a) performed a comprehensive assessment of different turbulence models for the simulation of the COPRA experiment. The authors concluded that the RANs models are adequate and recommended for the simulation of the natural convection behaviour of such internally heated melt pools. Moreover, DNS/LES are computationally expensive and time-consuming, which made it infeasible for simulations due to computational restrictions.
The k-omega turbulence model (two-equation turbulence model) is used in this study. Since the macroscopic properties of the prototypical material varies with the temperature, the three major conservation equations, namely, the mass, momentum, and energy are solved simultaneously (COUPLED Algorithm).

Reynolds-Averaged Navier-Stokes (RANS) turbulence models
In the RANS assumption, the parameters are segregated into the mean component (̄) and fluctuating component ( ′ ). The velocity, for instance, is decomposed as, With such an assumption, the mass and momentum conservation can be written as (Ansys ® , 2021), in which, is the velocity field, is the fluid density, is the Kronecker delta and is viscosity. A. Venkateshwaran et al. Additional equations for the solution of the Reynolds stress term (− ′ ′ ) must be included for the closure of the previous equations. The Reynolds-Stress turbulence models and the eddy viscosity turbulence models are usually utilised methods for solving this problem.

Eddy viscosity turbulence model and k-omega SST model
The Reynolds stresses could be approximated using a turbulent viscosity in eddy viscosity turbulence models. This can be stated as follows (Ansys ® , 2021): Hence, solving the turbulent viscosity is one of the most essential task in eddy viscosity turbulence models. Different turbulence models, such as Spalart-Allmaras, k-epsilon, and k-omega, can solve the turbulent viscosity based on different assumptions.
The turbulent viscosity in the high-Reynolds number is commonly estimated using the following equation in the SST k-omega turbulence model: in which, is the turbulent kinetic energy, is the specific dissipation rate, is the strain rate magnitude, * is the damping coefficient. The transport equations for the aforementioned turbulence model are as follows (Ansys ® , 2021), in which, is the generation of , is the effective diffusivity of , is the effective diffusivity of , is the dissipation of k due to turbulence, is the dissipation of due to turbulence, is a user-defined source term of , is the cross-diffusion term. The effective diffusivity for the − are as follows (Ansys ® , 2021), where and are the turbulent Prandtl numbers for and , respectively.

Solidification model in FLUENT
In this study, the solidification model in FLUENT is used to simulate the crust formation. The liquid-solid front is not explicitly modelled in FLUENT's solidification and melting model, rather, the porous zone method is used. In the porous zone, the porosity value is determined by the liquid percentage of the melt material. In the solidification model, the energy equation is as follows (Ansys ® , 2021): in which, where ℎ is the sensible heat, is the latent heat, ℎ ref is the reference enthalpy, ref is the reference temperature, and is the specific heat at constant pressure.
The following equations are used to calculate the liquid fraction (Ansys ® , 2021), where is the latent heat of the simulant material, The latent heat content can vary between zero (for a solid) and (for a liquid).
The mushy zone is defined as an area where the liquid fraction is between 0 and 1. Furthermore, if the liquid fraction is zero, the velocity will be zero. The mushy zone is considered as porous media in the momentum calculation, and the liquid fraction is assigned as the porosity for the porous medium. The force term,̄in the momentum equation can be written as . In such a mushy zone region, the porosity induced momentum sink is calculated as follows (Ansys ® , 2021): where is the liquid volume fraction, is a small number (0.001) to prevent division by zero, ℎ is the mushy zone constant, and ⃗ and ⃗ is the fluid velocity and solid velocity (caused by pulling of solidified material out of the domain), respectively.
The momentum sink terms that influences the turbulence model can be expressed as (Ansys ® , 2021): in which, is a small number to prevent division by zero, and is the turbulent quantities ( and ). The damping amplitude is measured using the mushy zone constant. The issues related to quick velocity drop during solidification and the unsteady convergence during the simulation could be caused by a higher mushy zone constant. Thus, a small mushy zone constant of 3000 is used.

Materials and methodology
In the current work, to simulate the COPRA test, the eutectic mixture of 50 mol% KNO 3 -50 mol% NaNO 3 was employed as the simulant material and the height of the corium pool was fixed as 1.9 m as per the experiment . The geometric parameter, boundary conditions and the material properties are explained in the following sub-sections: Fig. 2 shows the boundary conditions involved in this study. The dimensions are same as the vessel dimension used in the COPRA test facility but only two-dimensional geometry is considered here. For the sake of modelling simplicity, the focus is on the behaviour of the melt pool, hence only the melt pool is modelled for natural convection ignoring the air region above the melt pool. The isothermal boundary condition on the curved vessel wall acts as a surrogate for modelling the convective behaviour inside the vessel.

Simulation setup and boundary conditions
In the experiment, the portion above the melt pool is filled with air. In order to mimic the cooling effect of the air at the top, the boundary at the top is set as a radiation heat transfer which accounts for the heat loss between the air and pool in the top section of the domain. The emissivity for the corium pool is assumed to be 0.44 and the external radiation temperature is set as 300 K which is the ambient temperature of the experiment and same as the coolant water temperature. Adiabatic condition is assumed for the vertical wall of the pool.
The characteristic internal powers, which were calculated from the measurements along the curved wall (Luo et al., 2019), was about 47% total heating power for the eutectic salt, hence, 11.2 kW/m 3 was used as power density. The curved outer vessel wall is given an isothermal A. Venkateshwaran et al.  boundary of 300 K for the natural convection case and the interface of the vessel and the melt pool is thermally coupled. The temperature fluctuation along the curved wall outer surface could be maintained within ±1 • C thus creating the nearly isothermal boundary condition . The initial temperature of the simulant salt is set as 623 K (Luo et al., 2019), which is same as the pouring temperature of the salt in COPRA test. The fluid is also considered at rest initially.
Another set of simulations is performed, with a focus on the effect of forced cooling via the water flow. In the case of forced convection in eutectic salt, the input mass flow rate of coolant water is 8 m 3 ∕h  and the outlet is a pressure outlet boundary condition. Rather than using an isothermal boundary condition in the case of natural convection case, here, thermal coupling is implemented between the outer wall of the vessel and the coolant channel. The outer wall of the coolant channel is adiabatic with no heat transfer. The other boundary and initial conditions are same as the natural convection case as shown in Fig. 2. Table 1 shows the simulant material properties of 50% KNO 3 -50% NaNO 3 , that are referred from the prior studies (Luo et al., 2019). The material properties of the steel used in the vessel construction are also presented in Table 2. Higher order discretisation for better accuracy, and a mushy zone parameters of relatively low number which is explained in the previous section, i.e., 3000 was employed to get close results to the experiments.

Mesh independence test
A total of four mesh sizes were investigated for the mesh independence study and the simulation conditions (natural convection, no water flow) explained in 3.1 were employed. Fig. 4 shows the trend in temperature profile along a vertical line on the melt pool for different mesh sizes. The result converges when the total number of nodes in the mesh is greater than 75 000. The Table 1 Properties of simulant material (Luo et al., 2019).

Properties
Eutectic salt  percentage error between the two temperature curves converges from 1.95% to 0.40% at 75 000 nodes. Hence for all the cases, a mesh having 75 000 nodes (depicted in Fig. 3) is chosen. Table 3 shows the different geometric sizes implemented for this study along with the calculated internal heating density for each case. These three geometric cases are simulated in both natural and forced convective framework. The vessel used in COPRA experimental setup has a width of 0.2 m . In this study, two-dimensional geometry was chosen, hence the width of test section is not considered.  The volume of the melt pool is scaled proportionally to the radius of the vessel when varying the geometries. To calculate the internal heating density for the other geometries, the same ratio of heating power and cooling surface length was used. The internal Rayleigh number for each case is also calculated according to Eq. (19). The material properties of the salt at 350 • C are used for Ra ′ calculations. The same simulation methodology used for validation is extended for geometry with varying size, in order to understand its impact on the natural convection behaviour of the melt pool.

The effect of geometry
The relation between the internal or modified Rayleigh number, Nusselt number, and dimension of the geometry for an internally heated melt pool is given below (Lee et al., 2014): where, is the Nusselt number, is a numerical constant determined by the flow direction, ′ is the internal Rayleigh number, is the characteristic length (or height), is the radius of the vessel, and and are interpolation constants. The ratio of convective to conductive heat transfer in the melt pool is represented by Nusselt number ( ). The internal Rayleigh number, which is also referred to as the modified Rayleigh number, is pivotal in natural convective heat transfer flows that involve heat generation. It quantifies the internal heat source, thus the buoyancy effect (Rempe et al., 2005). This quantity is a modification of the existing Rayleigh number and is illustrated in Eq. (19).
where, is the Dammkohler number, is the conventional Rayleigh number, is the Grashof number, is the Prandtl number, is acceleration due to gravity (m/s 2 ), is coefficient of thermal expansion (1∕K), is the heat generation rate (W/m 3 ), is the characteristic length (m), is the kinematic viscosity (m 2 ∕s), is the thermal diffusivity (m 2 ∕s), and is the thermal conductivity (W/m K) (Rempe et al., 2005). The non-dimensional numbers are further given as: = .
The heat transfer via natural convection with internal heat generation is represented by the modified Rayleigh number( ′ ), that indicates the significance of buoyancy force in the melt pool. While Rayleigh number can describe the buoyancy and viscous forces in play, it does not account the internal heat sources. This necessitates the modification via the Dammkohler number and the internal Rayleigh number is of particular interest to understand the flow dynamics in such experiments (Rempe et al., 2005). This study evaluates this non-dimensional quantity as well for all the simulations performed.

Results and discussion
This section covers the results of the simulation of COPRA eutectic salt experiment using natural convection (without any water flow for cooling) and forced convection (with water flow for cooling). The comparison of results with the geometry of the melt pool is also discussed. The temperature contours, velocity profiles, formation of crust in the vessel, the time-series of turbulent structures in the melt pool and the Nusselt number-Rayleigh number correlation are presented. Fig. 5(a) and (b) portray the temperature field and crust formation, respectively, with the eutectic salt being used as the simulant material. The results of the present work were compared with the previously done study by Luo et al. (2019). Fig. 6 depicts the temperature varying along the height of the melt pool in this study and comparison with the experimental values (from COPRA ). The average percentage error between the experiment and numerical temperature profile is estimated to be 6.4%, whereas the maximum deviation is 9.15% These aberrations potentially arise due to the difference in geometrical consideration. While the original experiment was conducted in a 3D geometry (A quarter circular slice having a width of 20 cm), the numerical model is built upon a 2D geometry. Another plausible reason could be the internal decay heat. In the experiment, electrical heating rods are used at discrete zones of the melt pool, which are non-homogeneous, whereas the numerical model incorporates uniform decay heat across the pool using a source term. Moreover, the change in material properties with respect to temperature is a smooth function in the experiment, whereas piece-wise polynomials are implemented here. Other contributing factors could include numerical errors, inaccuracy in thermocouples, and so on.   Similar deviation from the experimental data can be observed in other numerical studies conducted by Wei and Chen (2020) and Luo et al. (2019). Overall, the present model reasonably captures various properties of the internally heated melt pool including thermal stratification, turbulence flow, and crust formation at the upper and lower region of the melt pool, respectively as recorded in the experiment .
Each relevant field (temperature, velocity, crust formation, and so on) is now discussed for the above simulation and more cases. All snapshots are taken when quasi-steady-state was reached in the simulations. For ease of perusal, these times are listed in Table 4.

Temperature field and temperature profile
The temperature distribution profiles for Case A (0.5 m pool radius), Case B (1.0 m pool radius), and Case C (2.2 m pool radius) via The entire melt pool is at 350 • C (or 623 K) initially. As time progress, the bottom region of the melt pool undergoes cooling due to the isothermal boundary condition on the curved vessel wall. Meanwhile, the top region of the melt pool is subjected to the radiation boundary condition, therefore, having a lower degree of cooling. As a result, a higher temperature is maintained in the top region. On  the other hand, a lower temperature is formed at the bottom and this creates thermal stratification. This stratification is further contributed to by the buoyancy effect of the melt pool. The hot eutectic salt floats at the top due to buoyancy whereas the cold molten salt settles down to the bottom of the pool and further cools to form a crust along the curved wall. This observation holds true for both natural and forced convection cases.
The maximum melt pool temperature for the smallest radius, Case A, is slightly lesser than the other two cases. With an increase in the geometric size of the melt pool, there is an overall increase in the thermal mass of the melt pool, that is the capacity of the melt pool to sustain high temperatures. This possibly results in a higher maximum melt pool temperature for larger melt pools. This can be attributed to the thermal inertia due to the size of the melt pools and this trend holds for both natural and forced convection cases.
Due to a noticeable temperature difference between the top region and bottom region accompanied by varying density of melt pool with temperature, strong turbulence and mixing in the upper region were observed. With the increase in the geometric size of the melt pool, an increase in turbulence is also observed. Fig. 9 shows the graph of temperature variation along the height that is normalised by the maximum pool height. The normalised height helps to maintain a same range (0 to 1) for all geometries, thereby comparison of temperature profiles for different geometries can easily be done. The profile in Fig. 9 can be divided into three regions. In Region I, up to a normalised height of around 0.15, the temperature sharply increases, which is attributable to the formation of crust and the cooling condition set at the curved vessel boundary.
The temperature rise in Region II is very gradual compared to the previous region. This can be explained by the fact that the velocity of melt pool in this region is observed to be very low. Subsequently, the molten salt temperature again surges in Region III, from a normalised height of 0.8, because of thermal stratification as seen in the temperature contour before. The temperature profiles follow a consistent pattern with the experimental data as shown in Fig. 9. Fig. 10 shows the similar trend and features of the forced convection. Region I is extended to a normalised height of 0.2 because the crust is thicker in all the cases of forced convection compared to natural convection. The deviation from the experimental data is expected and are explained previously in Section 4, it can also be observed that the Region I and III cover more portion of the melt pool in the forced convection model thus adding to the point that the thermal stratification is well established in the case of forced convection when compared to the natural convection furthermore substantiating the argument of considering forced convection model to provide better result than natural convection model. Fig. 11 shows the velocity profile for the Case A, Case B, and Case C depicting the influence of geometry on the motion of corium pool in natural convection. Fig. 12 depicts the same for forced convection.

Velocity and streamlines profile behaviour
The velocity profile follows similar trends of temperature as well. The maximum velocity tends to increase with an increase in the geometric size. High velocity is observed at the top of the melt pool because of the radiative boundary condition. Even though the top region of the melt pool has strong turbulence, the maximum velocity is seen along the curved wall especially above the crust region as the fluid rushes down the side of the curved vessel wall. This can be attributed to the cooling effect of the isothermal boundary condition. This is also evident in the vector plots presented in next section.
The temperature gradient creates the circulation of the melt pool. As this re-circulating flow encounters the crust region on the curved wall, it breaks down and creates vortices as observed in Fig. 11. These perturbations cause better mixing of the melt pool and create a more stable thermal stratification near the crust. This can be seen in Fig. 7 as the temperature distribution is uniform where these vortices are formed. The trend of temperature profiles in Fig. 9 is also largely stable in region two until it rises at around 0.8 normalised pool height. The velocity of these vortices increases with the increase in radius of melt pool. This directly affects the natural circulation intensity thus affecting the heat transfer characteristics along the curved wall. The maximum velocity at the frontier of the crust is found to be 0.002090 m/s, 0.002673 m/s, and 0.003993 m/s for Cases A, B, and C respectively. Similar magnitudes are also noted in the results of forced convection (Fig. 12).
In the case of forced convection (Fig. 12), turbulence is predominantly seen at the top region of the melt pool, which is also observed in COPRA-relevant prior studies (Luo et al., 2019. However, A. Venkateshwaran et al.     the same phenomenon is not observed in natural convection case as in Fig. 11. Thus, the simulation using forced convection provides considerably closer results to the literature and adequately represents the COPRA experiment. This is logical since COPRA does make use of a cooling medium, which is what is simulated in the forced convection case. The natural convection simulation, meanwhile, is indicative of a melt pool SAM strategy when water fails to flow and cool the vessel, therefore leaving the melt pool untouched. To better understand the relation between the velocity and crust formation, the overlap of liquid fraction contour and streamline for all the cases are presented in Fig. 13 for the natural convection case and in Fig. 14 for the forced convection case. It also provides an insight into the path of natural convection in the corium pool. The prominent feature in all the cases is the presence of vortical turbulent structures near the crust surface. This leads to the mixing of corium pool laterally. Similar findings, denoted as ''mixing region'', are also reported in Luo et al. (2019). The fluctuation in streamlines in the upper region of the melt pool is prominent in Case C, that is, the vessel with the largest radius. This can be also noticed in the velocity contours in Fig. 11 (see Fig. 14).

Evolution of eddies
The evolution of eddies formed in the eutectic melt pool is shown in Figs. 15, 16, and 17 for Case A, Case B, and Case C, respectively. A total A. Venkateshwaran et al.   Table 5) of velocity plots for each cases of geometry is exhibited. Due to the similarity in both the cases, only the forced convection case is shown for illustration purposes. The contours start from the moment of the peak temperature until the melt pool cools down to the quasi-steady-state presented.
It is evident from all three cases that there are two salient features of the velocity distribution that are consistent Firstly, is the formation of highly turbulent eddies at the top region of the melt pool. The evolution of these eddies can be structured into three stages.
1. The first stage of the evolution does not have any vortices at the top region. As specified initially, the velocity of the entire melt pool is initialised as zero at the start of the simulation. Thus, the first moment in Figs. 15, 16, and 17 do not show much pronounced fluctuations at the top region. As the transient simulation progresses, the circulation of fluids in the melt pool begins. This circulation is caused by the both gravitational and buoyancy effect in the fluid. 2. The next stage begins where the fluid gains momentum and circulates, the emergence of eddies at the top can be seen in Figs. 15 (snapshot 3 onwards), 16 (snapshot 2 onwards), and 17 (snapshot 2 onwards). Strong turbulence develops in all three cases, and it is more pronounced and larger in Case C compared to the rest.
A. Venkateshwaran et al. It is also evident from all the profiles that the eddies at the top region increase initially to a larger degree, and the intensity progressively reduces. This can be explained by the fact that as the temperature of the melt pool decreases the density and viscosity of the fluid increase. As a result, the fluid becomes more viscous, less chaotic, and turbulent. This causes more dissipation and the reduction of eddy size. This transition can be clearly seen in the images, such as the change in snapshots 3 to 4 in Fig. 16. This stage of velocity plots has the biggest coherent structures, whereas the subsequent moments show smaller eddies and diminish in time. 3. The final stage of turbulence at the top region of the melt pool is when the quasi-steady state is reached. Although the vortical structures oscillate and drift to different locations, the behaviour and size are persistent in time. The eddies at this final stage in Case C, Fig. 17, that is having the largest pool size, is more prominent compared to the rest. In Case A, Fig. 15, the eddies almost become negligible at the top. This is one of the geometric effects due to reduction in vessel size.
The second salient feature in these time-series plots is the downward flow along the cooling wall. A downstream flow along the crust's border starts at the cooling surface in the top half of the area. Since the cooling surface adjacent region has a comparatively strong cooling capacity, the temperature of the pool at the bottom is lower than that of the bulk, resulting in a higher density at the bottom. The downward flow of the melt pool with a relatively high velocity at the beginning will flow along the crust border, influenced by gravity. The flow is also slowed down by the shear force exerted by the crust surface and the viscosity of the fluid. The crust formed at the bottom hinders the downward flow. The evolution of this downward flow along with the formation of the crust is illustrated in all these profiles. Fig. 18 illustrates the vector plots at three stages in Case C, which corresponds to 1620 s, 6470 s, and 12 950 s time instances from Table 4. It depicts the direction as well as the magnitude of velocity vectors at evenly distributed points in the melt pool. As explained before, in the first stage, the eddies are disseminated and barely visible at top of the melt pool. These turbulent structures distinctly encompass the top of the melt pool in the next stage, however, it still has not reached the quasi-steady state, as a result, the size of the coherent turbulent structures progressively decreases. The final stage portrays the vector distribution at the quasi-steady state.
Furthermore, the peak velocity vector along the curved wall is at a low polar angle. It gradually ascends to higher polar angles and the magnitude decreases in the upcoming stages as the crust grows along the curved wall.

Crust distribution and heat flux
The liquid fraction contour is shown in Figs. 19 and 20 for natural and forced convection respectively. The crust thickness increases as the melt pool size decreases and reaches higher polar angles for a smaller size of vessel. This is true for both the natural and forced convection case. Another observation from Figs. 19 and 20, is the distribution of liquid fraction as the polar angle increases. The blue region denotes the crust (or solid phase) which steadily decreases as the polar angle increases for all the cases.
The formation of crust occurs along the curved wall as the melt pool is cooled by the relatively cooler isothermal curved vessel wall. Variations of crust thickness with respect to the polar angle for different radii are depicted in Figs. 21 and 22. It can be observed that the crust is thicker at bottom of the melt pool and gradually decreases along the curved wall as the polar angle increases. This trend is observed for all the geometries. As the vessel radius decreases the percentage of area occupied by the crust to the total area of the melt pool increases. This can be seen in Fig. 21, Case A (with a pool radius of 0.5 m), has the highest crust thickness ratio to the melt pool with more than 10% of the melt pool height being occupied by the crust, whereas, Case C peaks at less than 8%. This implies that the crust occupies larger area of the pool when the vessel radius is smaller.
In the case of forced convection in Fig. 23, Case A again has the thickest crust compared to the other two cases. The peak reaches a value of 20% of the melt pool which is twice the value of the natural convection discussed previously. Furthermore, Case C has its peak at 10% of the pool. Hence, the crust associated with forced convection is much thicker and more aligned to the experimental observations  as compared to the natural convection simulations. It can be noticed that the numerical model deviates from the experimental study. Similar observations are also recorded in previous studies by Wei and Chen (2020) and Luo et al. (2019), this arises mainly due to the consideration of natural convection as the predominant mode of heat transfer. This deviation is completely removed when forced convection is implied, thus bolstering the argument that the forced convection model can predict crust formation accurately.
The crust formation is under-predicted for Case C (R = 2.2 m) when compared to the experimental data from Zhou et al. (2018). Due to the shrinkage of salt during sudden solidification, a small air gap can be formed which has a thickness of 0.3 mm (Wei and Chen, 2020). The effect of the gap is not simulated for simplification.
Moreover, from Fig. 22, as the radius of the vessel decreases, the crust formation reaches higher polar angles. For Case C (vessel radius 2.2 m), which has the largest radius, the crust thickness reaches zero at a polar angle of about 70 • . Meanwhile, the crust thickness for Case A (vessel radius 0.5 m) is sustained at this location and reaches zero only at about 85 • polar angle. The same pattern is consistent with the forced convection as well and is shown in Fig. 24.
This can have a significant impact if this design was used in a smaller RPV (such as in a small modular reactor) and the water flow is interrupted. Although the chances of this occurring are small, we recommend that this should be included in the probabilistic safety assessment (PSA). Further research and experiments are also recommended for this.
The trend of normalised crust thickness in Case C, indicated by the red line, fits the experimental data better in the case of forced convection (shown in Fig. 24) than the natural convection (shown in Fig. 23). The red curve in forced convection commences with a sharp fall at a polar angle for about 7 • of polar angle followed by mild increase in thickness. This trend is also observed in the trend of experimental data points (see: Fig. 24).
The formation of crust along the curved wall behaves as a thermal insulator between the vessel and melt pool, thus impacting the heat flux. Consequently, the heat flux is reduced at low polar angles compared to the top region of the melt pool. As a result, the heat transfer rate is lower at the bottom region of the curved wall compared to the upper region where the crust is not present.
Another argument to augment the explanation of the heat transfer rates could be that the temperature difference between the bottom region and the wall is much lesser than the temperature difference between the top region and the wall, thus supporting the fact that the A. Venkateshwaran et al.          heat transfer rates are more near the upper region while it is lower in the bottom region.
In the experiments, the heat flux along the curved wall is calculated using (Zhang et al., 2016b), where, local is the local heat flux, wall is the thermal conductivity at the wall, is the temperature of the vessel wall, is the thickness of the vessel, is the vessel radius, and the subscripts and represent the inner and outer surface of the curved wall respectively.
In this study, the average heat flux of the melt pool is calculated using the surface area-weighted average, identical to the method followed in the experiment (Zhang et al., 2016b). The weighted average of heat flux over the entire curved wall is derived as: where, mean is the mean heat flux, local is the local heat flux from Eq. (23) and is the lateral surface area of the heating zone. In the experiment, the melt pool was vertically divided into 10 heating zones, each measuring 190 mm in height and the temperature is measured using the external and internal thermocouples installed on the vessel at these zones. These temperatures were used to determine the local heat flux at individual heating zones according to Eq. (23). Subsequently, the average heat flux, encompassing all zones, is calculated using the surface area associated with these zone according to Eq. (24).
The local heat flux along the curved wall is illustrated in Fig. 25 for natural convection and in Fig. 26 for forced convection. The heat flux in all the cases is often minimal and consistent in the region where the polar angle is less than 60 • . This is due to the thermal resistance provided by the crust between the melt pool and the vessel. Thus minimal heat transfer takes place up to 60 • . The results for the forced convection are more in line with the experimental data, especially for Case C as it has the same dimensions.
Reiterating from the previous discussion, the normalised crust thickness decreases to zero at a polar angle of 70 • for Case C. As a result, the thermal resistance is diminished beyond 70 • , leading to a steep increase in the heat flux, which can be seen for Case C in Fig. 25. The same phenomenon is observed for all three cases in forced convection. It can be astutely observed that the maximum heat flux increases as the size of the pool decreases.
The heat flux is minimal until a polar angle of 50 • prior to a surge in value to 20 kW/m 2 . The figure hovers at around the same value till 70 • followed by a steady increase to 32.5 kW/m 2 .

Nusselt number and Rayleigh number
The Nusselt number can be calculated for the upper and curved surface (downward) of the melt pool. Since the heat transfer at the upper surface via radiation is considerably less, for illustration purposes, Nusselt number is determined only for the curved surface denoted by ''downward'' Nusselt number( ) which is presented as (Zhang et al., 2016b), where, int is the temperature on the inner surface of the curved wall, max is the maximum melt pool temperature, is the pool height, is the thermal conductivity, and q mean is the mean heat flux from Eq. (24).
The internal Rayleigh number against the downward Nusselt number for three cases is shown in Fig. 27. The modified Rayleigh numbers ranges from 1.270 × 10 13 to 5.860 × 10 15 . For the geometry with a larger radius, the downward Nusselt number obtained from natural convection tests was higher. The Nusselt number and Rayleigh number derived from the experimental results of Zhou et al. (2018) are also shown in Fig. 27, and the respective correlation between these two dimensionless numbers is expressed as: = 0.0449 ′0.2556 (26) The internal Rayleigh number calculated according to Eq. (19) is the same for natural and forced convection. However, the Nusselt number varies for the natural and forced convection based on the heat flux according to Eq. (25).
The Nusselt number obtained in this present work agrees well with the correlation function from the experiment . Generally, the natural convection case is slightly larger than the experimental data, while the forced convection agrees better with the experiment's trajectory.

Conclusion
A CFD study on the effect of geometry size on natural convection characteristics of eutectic (50 mol% NaNO 3 -50 mol% KNO 3 ) melt pool was performed using ANSYS FLUENT 20.2. The simulation was performed using existing experimental and numerical data. The melt pool temperature, velocity contour, crust thickness distribution along the vessel, and dimensionless numbers (Nusselt and modified Rayleigh number) for various sizes of the vessel were investigated.
This numerical study was divided into two cases, viz., cooling via natural convection and forced convection. In the case of natural convection, the isothermal boundary condition was used rather than   The yellow markers imply the experimental data points from Zhou et al. (2018), while the black line is the correlation function from Eq. (26). tested. The data points in forced convection were slightly more proximal with the correlation function compared to natural convection. 6. The prominent differences between forced convection and natural convection were in crust formation, heat flux, and velocity contour. As for crust thickness, forced convection prevailed with a thicker crust and closer distribution to the experimental results. As a consequence, the heat flux profile also procured a better trend over polar angle compared to natural convection. Finally, the eddies are captured accurately in the case of forced convection, especially at the top region of the melt pool. Overall, forced convection has proven to be a better candidate to simulate the melt pool behaviour compared to natural convection.
This study provided a better understanding of thermal-hydraulic behaviour in eutectic melt pools of different geometry during postulated severe accidents.
Melt pool natural convection behaviour with two different simulant materials is complex. In the future, we recommend examining the influence of the shape of the vessel, material properties, and the impact of the coolant on the heat transfer characteristics of the melt pool.
The results of the study are in line with the historical way in which IVMR studies have been performed. The Loviisa-440 reactor was backfitted with the IVMR strategy for SAM system, and was proven to be safe. As reactors grow larger, the reactor pressure vessel sizes grow, and so do the uncertainties that come with it. Recent experiments have addressed this concern successfully as well.
This study provides a perspective on the other end of it, taking the data from the larger RPV-based experiment and looking at what would happen as the RPV size decreases. As the promise of SMRs increases, we believe that this study enhances the confidence in the SAM strategies based around smaller reactors.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.