Progress in modeling solidification in molten salt coolants

Molten salts have been proposed as heat carrier media in the nuclear and concentrating solar power plants. Due to their high melting temperature, solidification of the salts is expected to occur during routine and accidental scenarios. Furthermore, passive safety systems based on the solidification of these salts are being studied. The following article presents new developments in the modeling of eutectic molten salts by means of a multiphase, multicomponent, phase-field model. Besides, an application of this methodology for the eutectic solidification process of the ternary system LiF–KF–NaF is presented. The model predictions are compared with a newly developed semi-analytical solution for directional eutectic solidification at stable growth rate. A good qualitative agreement is obtained between the two approaches. The results obtained with the phase-field model are then used for calculating the homogenized properties of the solid phase distribution. These properties can then be included in a mixture macroscale model, more suitable for industrial applications.


Introduction
Molten salts coolants have been recently proposed as heat carrying media in nuclear [1] and solar concentration [2] energy applications. This is mainly because of their high heat capacity, low vapor pressure and chemical stability in high temperature and radiation environments [3]. In particular, molten salt reactors (MSR) [4] are a central subject of research in the nuclear --( -). The molten salt flows into the reactor cavity from the bottom at about 700°C after passing through the heat exchangers (see figure 1). In the core cavity of the reactor, 3 GW of thermal power are produced due to nuclear reactions. Consequently, the temperature of the molten salt rises approximately by 100°C. The molten salt exits the reactor cavity from the top and is then propelled by the pumps back into the heat exchangers. In case of a loss of cooling accident, the fuel salt temperature may increase beyond the acceptable limits due to the nuclear decay heat. In order to avoid causing damage to the reactor structures, a passive safety system has been designed. This passive system consists in a set of cold plugs located in the lower part of the reactor, between the heat exchangers and the core cavity. The cold plugs are made of solidified salt. The working principle is based on the melting of the cold plugs as the reactor temperature rises, which allows the fuel salt to drain into dedicated storage tanks. Since this is the main safety system of the reactor, the cold plugs have to be carefully designed, to substantially reduce uncertainties. For example, a typical design requirement is that the cold plugs should not melt during routine operation. On the contrary, they should quickly melt in case of a large temperature rise in the reactor.
The SWATH experiment is currently being developed at the LPSC laboratory (Grenoble -France) as part of the SAMOFAR project [6]. One of the objectives of this project is to theoretically and experimentally study the solidification phenomena of a molten salt in order to improve the design of the safety systems [7]. The molten salt selected for the SWATH experiment is a eutectic FLiNaK (LiF-KF-NaF, 46.5% 42.0% 11.5% --mol) salt. This salt does not contain nuclear fissile material. However, at the operation temperatures in the SWATH experiment, the salt has similar dimensionless numbers (Prandtl (Pr), Weber (We) and Lewis (Le)) than the fuel salt of the MSFR. This allows verifying the performance of the fuel salt solidification models with a ternary FLiNaK system in the SWATH experiments.
Due to the importance in industrial applications of eutectic alloys, various analytical and computational models have been developed for studying eutectic solidification [8,9]. As base theory, the well-known formulation of Jackson and Hunt [10] describes steady-state eutectic growth and relates the lamellar spacing of phases to the undercooling of the solid-liquid interface, for a known velocity of this interface. According to this theory, a simultaneous growth of the lamellae in binary alloys occurs at an optimum lamellar spacing. This optimum spacing corresponds to the minimum undercooling of the solid-liquid interface [11]. However, recent numerical studies for binary eutectic systems [12] and for ternary eutectic systems [13] have shown that the patterns formed in this lamellae may depend also on the solidification path. In particular, it was demonstrated that there might be several solidification structures that exhibit a stable growth with similar undercooling.
The computational modeling of eutectic solidification is challenging due to the presence of a liquid phase and two or more solid phases (depending on the thermodynamic phase equilibria of the system). Furthermore, the challenge of the numerical simulation of eutectic solidification is to obtain an accurate resolution of the dynamics of the interfaces formed by the junctions between the different phases. Table 1 presents a brief summary of some of the different methods used in the modeling of eutectic solidification.
The enthalpy method is one of the most popular approaches for numerical simulation of solidification [21]. This method consists in solving the energy conservation via the accumulated enthalpy. Then, the value of this accumulated enthalpy for each computational cell in the mesh can be used to track the interfaces. This method is interesting since it can be easily coupled-by operation splitting techniques-to the equations of conservation of species and momentum in the system [7]. Furthermore, it can be applied to large scales [22], like the one of the cold plugs. However, the solidification process resolved by this model is strongly dependent on the thermo-physical properties of the solidifying media, namely the thermal conductivity (which may be highly anisotropic in the solid phase), the solute diffusion coefficients and the spacing between the solid corps of the system in the mushy-zone.  [14] Macroscopic enthalpy method for eutectic solidification of binary systems. [15] Quantitative multiphase multicomponent phase field model for ternary alloy solidification. [16] Quantitative multiscale phase field model for large-scale ternary eutectic solidification. [17] Multiscale diffuse interface model solving the large scale by an averaged macroscopic system and the small scale by a phase field method. [18] Phase field method coupled with atomistic simulation models such as molecular and Monte Carlo methods. [19] Level set method for eutectic solidification of multicomponent systems. [20] Phase enthalpy method for the solidification of a binary alloy.
Furthermore, enthalpy models do not contain information about the distribution of the solid phases during the solidification process. These limitations will restrict the usage of this kind of models for the reliable design of the cold plug. In order to overcome these limitations, the solidification process needs to be studied with a mesoscale solidification model. This type of model allows calculating the thermal conductivity, solute diffusivity and corp-spacing values needed by the macroscale enthalpy model. This is the objective of the present work. The phase field method has been largely used for the study of the microstructure formation during solidification processes [23]. This method introduces a phase-field variable, e.g. f, function of position and time, to describe whether a computational cell is solid or liquid. In the case of a multiphase system, the phase variable can be also used to track the phases distribution in the solid (e.g. α solid, β solid, etc) [15]. The particularity of this method is that the phase variable allows to represent the various interfaces (e.g. solid-liquid interfaces and α solid-β solid interfaces) in a diffuse manner, so that no interface tracking is needed. For these reasons, the technique of phase field modeling will be adopted in this work for the study of the solidification process of the ternary eutectic system FLiNaK.
Recently, the solidification of the ternary metallic Ag-Al-Cu system has been studied by means of large scale phase-field simulations [24,25]. In these works, different brick-like patterns have been observed during the solidification process, giving rise to a spiral-growth in three dimensions. The obtained microstructure was analyzed and compared to experimental microphotographies with a good agreement. As noted by [26], the morphology of this microstructure is highly sensible to surface tension values. The surfaces tensions between phases in the FLiNaK system do not exhibit large differences (see 3). Therefore, the brick-like structure and the spiral-growth are not expected in this system. Furthermore, the solid phases are expected to grow parallel to the temperature gradient (i.e. without a tilt angle) as observed in an illustrative case for the similar LiF-NaF system in figure 2. Furthermore, the influence of the solidification path in the solid structures formed has been studied in [13,27]. It has been observed that the morphology of the solid structure will be sensible to changes in the volume fraction of the solid phases. For the solidification conditions of the MSFR, there are no apparent reasons to expect changes in the volume fractions of the phases and, hence, this phenomenon is gloss over in the present work.
Other ternary systems have also been experimentally studied, such as the metallic-alkali ternary system Mg-Al-Cu, which presents an intermetallic compound formed during solidification, determining the development of a regular/irregular eutectic structure. In particular, it has been observed in experiments [28] that the changes in the solidification path affect the regularity and twinning of the phases formed, rather than the distribution of the phases themselves. For the FLiNaK case, this effect is negligible, since there is a perfect partition of the phases during solidification. Very few works are available in the open literature on the development of numerical models for the solidification of molten salts. For example, the solidification of a molten salt in a pipe has been studied by means of numerical simulation using a macroscale enthalpy formulation in the work done by [29]. However, this macroscale model makes some important approximations. In the first instance, it considers the thermal conductivity in the formed solid as being isotropic and uniform, which may lead to an inaccurate temperature field in the system during the solidification process. In addition, the arm spacing of the solid corps in the mushy-zone is not considered, which make it difficult to estimate the momentum sink terms These two approximations can have a substantial impact on the model results and, therefore, ought to be improved. The aim of the present work is to address this research gap. In fact, this article aims at improving the understanding on the development of the microstructure during the solidification of a molten salt, using afterwards this information to increase the reliability of macroscale models.
This article is divided into two parts. In the first part, a semi-analytical model is developed for calculating the microstructure for the eutectic FLiNaK system at a steady solidification rate. The structure is computed by minimizing the undercooling of the solidification interface. The model equations are solved with a semi-analytical method, for guaranteeing that the predicted solid microstructure results from the model assumptions and not from the implemented numerical method. The implications on the solidification dynamics of the calculated solidified structure are then discussed. In the second part of this work, the solidification of eutectic FLiNaK is analyzed by means of a phase field model. Further insights in the solidification dynamics are identified and analyzed. Finally, some of the solid macroscopic properties such as thermal conductivity, solute diffusivity and corp-spacing parameters are calculated. These properties could then be used as inputs in a macroscale enthalpy model.

Description of the system
Molten salts used for industrial purposes are usually eutectic systems. These systems allow lowering the melting point, reducing the risks of unexpected solidification and improving the operating margins [30]. Little experimental data exists on the study of solidification of molten salts. However, the eutectic solidification process has been studied in other systems such as metallic [31], organic [32], ceramic [33] and glass [9] systems, amongst others. Therefore, part of the eutectic solidification-modeling techniques developed for these mixtures can be used for molten salts. Nevertheless, the particularities of the solidification of molten salts have to be taken into account when adapting the models developed for the aforementioned systems.
In order to decrease the vapor pressure in the systems operating with a molten salt at high temperature, salts composed by ionic bonds between halogens and alkali metals (e.g. LiF, KF, NaCl, etc) are chosen [34]. Consequently, the solid-solid and solid-liquid interface tensions in these systems are high » [38], which explains the non-faceted and regular structures normally observed during the solidification process of molten salts. For illustration purposes, figure 2 shows the solid structure obtained for an eutectic mixture of LiF-NaF (61% 39% mol) in an experiment performed at LPSC-Grenoble.
In the present section a theoretical study for the solidification phenomena of the ternary eutectic FLiNaK (LiF-KF-NaF, 46.5% 42.0% 11.5% --mol) system is presented. The ternary phase diagram for this system is shown on the Gibbs simplex presented in figure 3 [39]. The solidification temperature for the eutectic composition, known as the eutectic temperature, is 727 K. During the eutectic solidification process, three phases are produced LiF, KF and NaF, referred from now on as α, β and γ respectively. As observed in the phase diagram, no solubility is expected between phases. Therefore, during the solidification process, a perfect partition of the molten salt in the three solid phases is obtained. While the solidification process will more often occur under external convection conditions (for example on the surface of the plates of the heat exchangers of the MSFR), for simplicity, only heat exchanges due to thermal conduction are considered near the solid-liquid interface. The modeling efforts in the present work are focused on the microscopic structure modeling and not on the influence of the external convection rate on the solidification process, which will be the subject of future developments. Finally, in our applications the volume of the hypothetical solidified fuel salt will be almost negligible in comparison with the total volume of the molten fuel existing in the MSFR. Hence, a permanent renewal of the liquid in front of the interface by the external free stream convection can be assumed. Therefore, the morphological changes in the solid phase microstructure due to modifications of the volume fractions of the phases caused by a non-constant composition of the liquid in the solidification front (solidification path effects [13,27]) are neglected.
Starting by the ground theory of Jackson and Hunt, the solidification of binary alloys have been widely studied in the literature [10]. However, the dynamics of the solidification process for a binary alloys are different to that of a ternary alloys [40]. Therefore, a different theoretical treatment is needed. Consequently, in the following section, a theoretical framework is proposed for the study of the solidification process of eutectic FLiNaK. Furthermore, this model is then used for predicting the shape of the solid structures formed during the solidification process by minimizing the undercooling in the solid-liquid interface. Ternary phase diagram for the KF-LiF--NaF system [39]. The eutectic point is at a composition of (46.5%-42.0%-11.5%)mol and the eutectic temperature is at 454°C (727 K).

Methodology
The initial system is assumed to be molten FLiNaK, in a 3D domain (V R 3 Ì ) that is assumed to be infinite in the x y , directions and semi-infinite in the z direction, which is the growth direction for the solid (see figure 4). Initially, the system is in liquid state, at a uniform eutectic composition over the whole domain and at the eutectic temperature. At the initial time, heat is extracted at a constant rate from the inferior part of the domain (z = 0) and the solid starts to grow in the z direction.
Initially, the solid will grow erratically in the z-direction, due to the non-uniformity caused by the nucleating and coarsening of the initial solid nuclei formed and the initial solid structure rearrangements [41]. As the solid grows, the heat extraction process becomes less efficient because the solid imposes an additional thermal resistance between the heat sink (located at z = 0) and the solid-liquid interface, from where heat is being extracted. Moreover, the solid phase thermal resistance is larger for the (x, y) positions in which the solid have grown further along the z direction due to the initial perturbations. Consequently, after the solidification front has grown a certain distance from the initial plane, the initial perturbations are dumped and the solid-liquid interface continues to grow in the z-direction as a planar front [42]. This is the initial condition from which the model is built.
As the solidification process is taking place, a redistribution of solute occurs in the liquid due to the difference in the composition between the liquid and the solid phases [43]. For example, as the α-phase is solidifying, NaF and KF are rejected by the solid phase, while LiF is absorbed by it. Therefore, for the energetic cost of solute migration, the configuration of less energy in the solid phase is the one in which the three phases of the system shown in figure 4 are infinitely close to each other, i.e. having a negligible thickness. However, due to crystallographic differences of the three solid phases, there is an energetic cost induced in solidifying the three phases next to each other. The energy for the formation of interfaces between two phases n n¢ in the system can be computed as the specific energy cost of forming an interface between two phases n n¢ , multiplied by the total area of the interface (A nn¢ m 2 [ ]). Summing over all types of interfaces in the system, the total energy cost associated to the interfaces in the system is Therefore, for reducing the energetic cost of solidification, the system will tend to minimize the interfacial area, thus maximizing the thickness of the formed solid phases. The competition between these two phenomena explains the formation of regular structures of finite length in the eutectic solidification process [44]. Therefore, in this section we study an infinite system assuming periodicity in the structures formed during the solidification process and using a unit cell to represent the smallest repetitive pattern of this infinite system (see figure 4).
The previously discussed energetic costs of solute migration in the development of the interfaces, is traduced in the system as an undercooling of the solid-liquid interface below the eutectic temperature. For a given phase ν ( , , n a b g = ) the defect in the interface temperature with respect to the eutectic temperature (T E ) can be calculated by the Gibbs-Thomson relation [45], formulated for the present system as, where c X is the concentration in the liquid adjacent to the interface of component X (being is the liquidus slope evaluated at the eutectic composition. Furthermore, T L l E g G = n n n is the Gibbs-Thompson coefficient, where l g n is the solid-liquid surface tension and L ν is the latent heat of fusion per unit volume for the ν phase. The relevant parameters used for the present study are given in table 2. The concentration field (c X ) in the liquid is described by a stationary diffusion equation resulting from the conservation of mass of the component X in the liquid. The description of the system is carried out on an inertial reference frame moving at the speed of the solidification interface (v) in the z direction. Assuming for simplicity a constant diffusion coefficient for the three phases (D), the equation can be written as follows: where D X is the diffusion coefficient for the component X. For the studied unit cell domain, shown in figure 4, periodic boundary conditions are assumed in the x-y directions, c y z c L y z 0, , , , , 3 For the z direction, it is assumed that as z  ¥, the concentration goes to the original eutectic concentration c c X X E  . For z=0, following the Jackson-Hunt theory, the solidliquid interface flux condition is given by the Stefan's boundary condition [53] approximating the composition of component X in the liquid to the eutectic composition, , and defining the dimensionless parameters , the dimensionless problem for the diffusion of component X can be written as follows: The conservation equation ( where X nm are the coefficients of the Fourier expansion that have to be determined. By inserting the proposed solution (11) in the general equation (6), the decay coefficients q nm are obtained as follows: To determine the coefficients X nm of the Fourier series, the proposed solution (11) is inserted in the Stefan's boundary condition (10) yielding, Next, the dot product with norm is defined as, Taking the scalar product on both sides of equation (13) gives, Working this expression, noting that the indexes (o, p) are dummy ones that can be replaced by another set of dummy indexes (m, n) and that A c = 1, the Fourier coefficients can be determined as follows: Finally, as displayed in figure 4, the unit cell area is composed independently by the three solid phases α, β and γ. Therefore, the integral over the area in the right-hand side of equation (15), can be split over the area of each of the phases, where each individual phase contains a component X, resulting in, Once the distribution of the phases in the unit cell is known, the Fourier coefficients are uniquely determined. If experimental information of the distributions of phases in the solidified structure is available, equation (16) can be used for determining the concentration field in the liquid in front of the interface. However, to the authors knowledge, this experimental information is not yet available for the FLiNaK eutectic system. Therefore, a different approach is taken. The average undercooling of the interface is calculated for the unit cell and, by minimizing this undercooling, the solid structure is computed.
As expressed by the Gibbs-Thompson relation (1) [45], the undercooling of the solidliquid interface can be expressed as the sum of the constitutional and curvature undercooling.
, the Gibbs-Thompson relationship can be expressed in dimensionless form as, ). The average constitutional undercooling of the solid liquid interface can be obtained by the weighted average of the constitutional undercooling over each of the phases. This can be written as follows: where A n is the area occupied by the ν-phase in the unit cell. In order to compute the curvature undercooling, the curvature of the solid-liquid interfaces needs to be determined for each of the phases in the system. The curvature for each component is defined as the divergence of the normal unit vector to the solidification interface [54]. In the present case, following this definition, the curvature for a phase ν containing a component X is proposed to be computed as . Developing this vector relation in the x-y directions, the curvature can be computed by This curvature is numerically calculated on the x-y plane for each of the components X and its corresponding phase ν. The mean curvature undercooling of the interface is then determined as, Finally, the general undercooling of the interface can be computed as the sum between the constitutional and curvature undercooling, The solid structure is determined as the one minimizing the undercooling of the interface. The processes of optimization is carried out via a feasible direction method [56], in which the objective function to be minimized during the optimization method is the interface undercooling as a function of the spatial distribution of the solid phases, The constrain of the conservation of species is imposed in the concentration field is imposed as The overall optimization process is displayed in figure 5. Initially, a random distribution of the solid phases, satisfying the constraint in the concentration field (23), is proposed. With this distribution of the solid phase, the Fourier coefficients for the concentration field in the liquid (16) are computed, thus allowing obtaining the concentration field. Then, the undercooling of the solid-liquid interface is computed (21). In particular, the curvature undercooling (19) is numerically computed using a fourth order ENO-Pade scheme [55] for reducing potential integration errors. Afterwards, a feasible direction of optimization for each of the phases is computed [56] and the solid phase distribution is actualized to the new calculated value. This new solid phase distribution is used for recalculating the Fourier coefficients (16) and the iterative process is carried on. The area weighted difference between the iterations n 1 + and n in the solid fields is defined as, The iterations are carried out until 10 5  < -.

Results and discussion
The iterative procedure was carried out using an uniform mesh discretization of 1000, 1000, 5000 The results obtained for the solid phases distribution are shown in figure 6. Two plates of KF and LiF are predicted in the sides of the domain, whereas the NaF phase is predicted to grow in a column shape between the other two phases with a lens-like cross sectional area. Considering the periodic boundary conditions imposed in the x-y directions for the unit cell, the prediction for the formed solid field results in long plates of KF and LiF of similar thickness, which grow in the z direction. The center of these plates are in the sides of the unit cell. For the NaF phase, the growth is predicted along the z direction and between the plates of the LiF and LiF phases. The lens shape cross section x-y plane corresponds to the structure that minimizes the surface energy of the phase respecting the Young-Dupre equilibrium at its boundary.
From the solute migration point of view, one would expected that the phase with less concentration would be located between the other two phases. This is because the migration process KF and LiF from the γ solid phase will consume less energy if the migration distance is as short as possible for both solutes. The reason why the LiF and KF phases solidify as long plates is a results of the large surface energy involved in the formation of the interfaces between these two phases. Therefore, the solidification process tends to minimize the surface between the LiF and KF phases, arriving to a unique straight plane perpendicular to the x-y plane.
The concentration field for each component in the liquid adjacent to the solid-liquid interface c x y , , 0 X * * ( )in shown in figure 7. When a phase ν is solidifying, the component X, associated to that phase ν, decreases its concentration next to the solid phase and increases it next to the other two phases. This is because the component is absorbed by the phase ν and rejected by the other two. For example, we observe for the KF concentration in the liquid that the concentration of the component decreases close to the β solid phase and increases close to the other two. In addition, the minimum concentration of the KF component is found at the left boundary of the unit cell, corresponding to the center of the KF plate. A similar behavior is observed for the other two components. There are two important facts resulting from the analysis of the concentration fields. First, it is observed that concentration in the liquid next to the interface is smeared out by the diffusion process in the liquid, resulting in a maximum variation of the concentration of the components next to the solid-liquid interface to be of about 5%  . This fact is important, since it will limit the ulterior design of safety devices based in the rapid solidification of the molten salt. Secondly, it is observed that the variations in the concentrations of the KF and  LiF components in the liquid are almost independent of the third NaF component present with less concentration at the eutectic point. Therefore, the grow dynamics of the solid system can be described with relatively good precision by considering only the KF and LiF phases and neglecting the third one. This is a valuable statement from the system design point of view, since the solidification of a binary alloy is, in many aspects, simpler to treat than the solidification of a ternary alloy. The cross-diffusion process occurring in the liquid, starting from the release of solute in two of the phases and ending with the absorption in third phase is shown in figure 8. For the KF, the maximum concentration is achieved in the middle of the plates of LiF phase and the minimum in the middle of the plates of the KF phase. Therefore, the average cross-diffusion migration process of the KF component occurs from the center of the LiF phase to the center of the KF phase. In addition, due to the velocity imposed on the z growth direction coming from the evolution of the solid-liquid interface, the iso-concentration profiles are observed to be advected in the growth direction. The exact reciprocal process happens for the LiF component, in which the cross diffusion process take place from the center of the KF phase towards the center of the LiF phase. For the NaF component, the cross diffusion process occurs from the corners of the unit cell, in which the concentration of NaF is maximum, towards the center of the cell in which the concentration is minimum.
As the two-dimensional concentration profiles are studied for larger distances from the interface (i.e. increasing along the z* dimension), it is observed that the concentration profiles are homogenized towards the eutectic composition of the components. In addition, it can be noted that the homogenization is faster for NaF than for the other two components. This is because the mean migration distance for this component is smaller than for the others. Nevertheless, for z 1 * = , the relative defect of concentration c x y , , 1 for all (x, y). Therefore, for practical applications in the solidification of this system, it can be assumed that after z 1 * = the concentration field in the liquid is equal to the original eutectic concentration.
The undercooling for the unit cell interface is analyzed on figure 9. The constitutional undercooling is maximum for the right and left borders of the system, corresponding to the middle of the plates of LiF and NaF, respectively. This can be explained by the fact that these are the regions from which the released solute have to migrate the maximum distance to be absorbed (or inversely the absorbed solute have been expelled from further away). The curvature undercooling is higher where the concentration gradients are higher. Therefore, the curvature undercooling is predominant in the bulk NaF phase and next to the contact points of the LiF and KF phases, where large concentration gradients exist. Even though, the surface tensions for the system are high, the curvature undercooling is reduced to 1/2 of the constitutional undercooling for the calculated solid phases distribution.
Finally, the total undercooling is computed as the sum of the constitutional and curvature undercooling. It is observed in figure 9 that the total undercooling is approximately constant for the three phases. This is expectable since large temperature gradients in the interface will be smeared out by heat conduction parallel to the solid-liquid interface. The area averaged total undercooling, defined by equation (20), recasting the interface temperature to dimensional form, is −15.172 K.
Although the present analysis is useful to understand the growth dynamics of the solid phases of a eutectic molten FLiNaK system, there are some underlying assumptions that limit its applicability to real scenarios. First, the phases are assumed to steadily grow parallel to the z direction, which in most cases is not the actual situation. As observed, the solidification process for a ternary eutectics is a highly nonlinear process and, due to various instability mechanisms, the phases are continuously rearranged as the solid phases growth [40]. In addition, it have been assumed when deriving the interface boundary condition that the concentration in the liquid next to the interface is equal to that of the eutectic composition. While in the present model it was readily observed that the diffusion process in the liquid modifies the concentration profiles. As shown in figure 7, the concentrations of the various components in the liquid near to the solidification interface are still strongly dependent of this initial hypothesis. Furthermore, it has not been assumed the presence natural convection due to temperature and concentration gradients in the liquid. This may be a keystone effect in the dynamics of the growth process since natural convection could have a large impact on the concentration fields in front of the interface. Finally, the model describes the solidification field once the system has arrived to steady growth conditions, assuming the interface to be a planar front. However, it does not contain information about the initial stages of the solidification process. Therefore, for surmounting some of these issues, a more precise study of the FLiNaK system is necessary. This is performed via a phase-field method in the next section.

Methodology
The conserved quantities describing a solidification process are the concentration of species, the energy and the linear momentum in the system [43]. Therefore, the numerical solution of this set of equations is, in principle, sufficient for the description of the solidification process. However, from the numerical point of view, this type of simulations are strongly computational demanding, since they require an accurate tracking of the solid-liquid interface [57]. Therefore, during the past decades diffuse interface models have gained popularity [58] and, among them, the phase-field method is being widely used for the studies of the solidification phenomena [59]. In the present work, a thermodynamic consistent phase-field model is propose for studying the solidification process of the ternary FLiNaK system.
As previously discussed, the solidified 3D domain (V R 3 Ì ) can be divided in independent volumes V a , each containing one of the solid phases. A phase-field f n function, t V R x, : 0, 1 f´ n ( ) [ ], is introduced for describing the solid field ν at a position V x Î and at a given time t R Î . When the value of the phase-field t x, f n ( )is equal to 1 it means that the solid phase ν is present at the position x at time t and when it is equal to 0 that is not. The values of the phase-field between 0 and 1 are related to interfaces between different phases. Furthermore, a phase-field variable l f can be introduced for describing the liquid phase. A vector is constructed describing the phase-fields t x, , , , ). The constraint 1 is imposed on this vector, in order for this vector to represent a valid thermodynamic state everywhere in the system. In a similar form, a concentration vector is introduced for describing the concentration field of each of the components of the system t c c c c x, , , The equations for the evolution of the phase-field and concentration variable are derived from a grand-potential functional Ψ [60], following the formalism developed in [15], which is a function of the phase-fields F, the concentration c and the temperature T, where the first term in the rhs T c , , y F ( )is the grand-potential bulk density and the terms in the parenthesis are the gradient energy density a ,  F F  ( )and the potential energy density 1  w F ( ), ò is a length scale parameter associated to the interface width. The last two terms are introduced for representing the energetic cost associated to the formation of interfaces between the phases.
The driving force for solidification is the dynamic reduction of the free energy in the system. Using an ideal solution model [61], the free energy density f T c , , where L i is the latent heat of melting of component i, r n is the density of the phase assumed constant, T E is the eutectic temperature of the system, R is the universal constant of gases and h 3 2 ) is an infinitely continuously differentiable function used to describe the transition between the liquid to solid for the phase ν. Since in the FLiNaK system each of the phases is associated to an individual component, the melting heat for a phase is equal to the melting heat of a component and, therefore, L i is independent of the phase. Furthermore, the free energy of the liquid phase can be expressed as f T RTc c c, l n following the ideal solution model. The bulk grad-potential density can then be derived as the Legendre transformation of the free energy density.
However, as initially noted by [62] for binary alloys and [15] for multicomponent systems, if the grand-potential is dependent on the concentration field, the surface energy of the interfaces will be a function of the bulk free energy field in the system. This means that the energy stored in the interfaces will be dependent on the model used for the bulk free energy density and the thickness of the interface, leading to non-quantitative results. However, as noted by [15] , where ij d is the Kronecker's delta function and the vectors , the second order approximation for the free energy (26) can be expressed as following, The chemical potential can then be written in the following form, By using equation (28) an straight forward relation can be founded between the concentration and the chemical potential, which proves to be useful when deriving the equations for the chemical potential conservation. By replacing the concentration by equation (29) in the quadratic approximation of the free energy (27) and by defining the bulk grand-potential as the Legendre transformation of the free energy , the bulk grand-potential can be expressed in the following way, The grand potential of the system can then be recast using the chemical potential in the following form, The gradient energy density can be expressed as, where g nn¢ is the surface tension between phases n n¢ and q is a generalized gradient vector that is different from 0 at the solid-solid and solid-liquid interfaces. This energy is associated to the necessary extra energy introduced by the gradient of phase-fields in the system. The potential energy density of the interfaces can be expressed by a double and triple obstacle potential as, 16 . The first term in the rhs is associated to the energy accumulated on interfaces between two phases and the second one to the interfaces between three. This last term, allows avoiding spurious third phases on binary interfaces. For avoiding violations on the constraints of the phase-fields ( 1 during the evolution, the potential energy is taken to infinity w F = ¥ ( ) (a computationally large number) at the points where the constraints are not satisfied. This energy is associated to the extra energy necessary for forming a new double and triple interfaces caused by the surface energy between the phases of the system.
For describing the evolution of the phase-field variable a relaxed Allen-Cahn equation [63] is used, where t n is a time relaxation parameter introduced for reducing the numerical time in the solution of the system with respect to the real time of evolution and M f is a mobility parameter for the phase-field. Inserting equation (25) in (34), the following equation is obtained for the evolution of the phase-field f n , where λ is a Lagrange multiplier introduced in the equation in order to satisfy the constraint . In a similar form, an equation can be derived for the evolution of the concentration field based in the grand-potential. However, since the concentration needs to be conserved during the evolution, a Cahn-Hilliard equation [64] is used for the description of the evolution concentration field, where c t is a relaxation parameter for the time evolution of the concentration field and D Ci F n ( ) is the diffusion coefficient for the concentration field i in phase ν. Introducing equation (25) , is a Lagrange multiplier introduced analogous to the phase evolution case in order to satisfy the constraint c 1 on the concentration field. Since the decoupling between the interface energy and bulk thermodynamic potentials requires a formulation based in the chemical potential an evolution equation for the chemical potential has to be determined. By using the chain rule and equation (37) it can be noted that Rearranging the above equation and by using relation (29) for the concentration, the evolution equation for the chemical potential is derived in the following form, Usually an anti-trapping current is introduced in the chemical potential equation in order to account for an artificial trapping of solute in the solid by the diffuse interface [65]. However, in the present case, the coupling introduced by the anti-trapping current was observed to introduce artificial numerical oscillations, coming from the large surface tensions involved in the FLiNaK system. Therefore, it was not considered and the artificial numerical diffusion of the interfaces was systematically controlled for avoiding solute trapping effects to affect the system evolution.
Analogous to the previous cases, the evolution for the temperature field can be derived from the grand potential. This yield, where Cj k is the thermal diffusivity coefficient for phase ν. The general diffusion coefficient is defined as the average of the eigenvalues set of the concentration diffusion matrix D D ave , Ci l = n ( ). In a similar form, the average thermal diffusivity is defined as   T  T L  c  c  , , , , is the Lewis number. By non-dimensionalizing the equations of evolution for the phases (35), the chemical potentials (39) and the temperature (40), considering the previously introduced dimensionless variables, the following dimensionless system of equations describing the system evolution is obtained (note that forced or natural convection in the liquid has not been considered), The list of simulation parameters for the FLiNaK system are shown in table 3. The interface surface tension values have been indirectly extrapolated from experiments and may be subjected to errors. However, a sensitivity study has been performed by changing these values and the results, presented in the next section, were not observed change substantially. Therefore despite the uncertainties existing in the values of some of the physical parameters, the set of values provided in table 3 can be considered as a reasonable starting point for the study of the solidification dynamics of the FLiNaK system.
The domain and the boundary conditions used for the simulations are shown in figure 10. In the x and y direction periodic boundary conditions are assumed, the length of the domain is 10 non-dimensional units in both directions (as the solidification process have been observed to be diffusion limited this will allow to obtain approximately 10 phases in each direction). In the z direction, a uniform heat flux (extracting heat from the system) is assumed at z 0 * = while far field boundary conditions are used at z L z * * = . An initial gradient in the phase field variables is imposed at z=0 for initially perturbing the system. It is assumed no concentration gradient in the z direction at z=0 and, therefore, 0 Far field boundary conditions are as well used for the phase variables and the chemical potentials at z L * = . The numerical simulation is stopped before the far field boundary have a significant effect on the solidification process.
The domain is decomposed on a structured regular mesh of 1000 1000 10 000( ) elements in the x y z , , ( ) directions. The initial concentration is determined by a Voronoi tessellation method [69] following the eutectic composition for each of the phases. With the ll l l l l l l g g g g g g g g g g g g g g g g  generated concentration field a 2D phase-field simulation is performed in the base plane (z 0 * = ) [70] in order homogenize the phase distribution. This allows to reduce the sensitivity of the system to the initial conditions. In order to increase the speed in the simulations, a pseudo-spectral method is implemented, where the x-y are solved in the Fourier space and the z direction is solved by a compact 4th order finite difference scheme [71,72]. This last aspect, allows to reduce the numerical diffusion of the solid-liquid interface, decreasing the artificial solute trapping and limiting the error introduced by not using an anti-trapping current. Additionally, an adaptive mesh refinement box technique is implemented with an octree method [73], in order to increase the mesh resolution close to the solid-liquid interface. The simulations were done in 64(@2.5 GHz) cores and the total simulation time was of approximately 48 h.

Results
The results for the obtained solidification field are shown in figure 11. A cooperative growth is observed between the three solid phases. As previously predicted by the semi-analytical model, it is observed that, once for steady grow conditions, the α and β phase solidifies as long plates, whereas the NaF phase grows as oval-shaped columns between the other two phases. During the solidification process, fibers of the solid phases are formed in the direction of grow. Going forward in the z* direction, the fibers merge due to a reduction in the interface gradient energy or split due to a reduction in the bulk free energy (the driving force of solidification). The formed structure of fibers is observed to be continuous, i.e. there are no branches that suddenly stop growing or are suddenly created. In order to understand the dynamics of the growing process, it is proposed to study the process from the point of view of the individual phases.
The grow process for the LiF a  ( ) phase is shown in figure 12. Initially, when the equilibrium have been reached with the 2D phase model in the base plane, the α phase is distributed in small clusters with convex shape. This sort of shapes are formed because the diffusion process in the 2D plane is overestimated with respect to the real one occurring in three-dimensions, which have not been yet considered. The system searches, therefore, to minimize the driving force of solidification. This situation is similar to what occurs during the coarsening of solid nuclei next to a wall, since the diffusion processes are helped by the adsorption and desorption of solute in the walls and the presence of the wall limits the influence of the surface tension.
As the system grows in the z direction, the initial reduction in the driving force of solidification caused by the presence of small nuclei occludes initially the effect of the reduction in the gradient and potential. This is, for reducing the gradient and potential energy, the fibers of the α phase will tend to merge and coarse as they grow in the z direction. Nevertheless, this coarsening will increase the driving force of solidification. Since, initially, the driving force and the chemical potential are low for the system in equilibrium; the system has a big inertia for changing its configuration. However, as the solidification progress further in the z* direction, the chemical potential rises as the initial equilibrium distribution is perturbed. This results in the system actively minimizing its interfaces and going to the previously observed thin-plate structure.  While the system grows further in the z direction, it is observed that a perfect plate is not formed, as expected by the reduction in the interface energy. The phenomenon is originated by the presence of oscillatory instabilities in the growing process. While growing upstream, the α fiber may reduce its diameter, while reducing the driving force, due to high differences in the concentration fields in the liquid in front of the fiber with respect to the liquid that surrounds it. As the fiber thins, the sink of solute will be reduced with the size of the fiber, causing the external solute field to be more rapidly homogenized by diffusion processes. Therefore, afterwards, the fiber will coarse back. This type of oscillatory behavior is present over the whole domain of the FLiNaK system, since it should be continuously triggered by the high values in the interface tensions.
The analysis of the grow of the KF b  ( ) phase is similar to the one of the α phase. However, the growing process for the gamma phase is different, since the phase solidifies as columns between the plates of the two other phases. In addition, this columns have an ovalshapped cross sectional area similar to the one predicted in the semi-analytical model. This process is shown in figure 13. Initially, when solving the phase-field model in the 2D base plane, once again, the diffusivity of the solutes is overestimated with respect to the 3D case and a disperse phase with high connectivity and approximately equispaced clusters is formed on the base plane. As the solid starts growing in the z-direction, the columns start to bend to the space between the α and β plates, since these are the locations that reduce the driving force of solidification. Once the columns migrate to the space between the plates, they will form clusters, reducing the gradient and potential energy of the interfaces. As in the previous case, these clusters do not have a perfect lens cross-sectional area because of the presence of an oscillatory behavior in the system.
In order to extract condensed information from the grow process of the solid phases, nearest neighbors statistics are used [74] (see figure 14). It can be seen that the amount of nearest neighbors in the system is reduced for all the phases as solidification progresses in the z direction, to finally reach a constant value in steady grow conditions. For the α and β phases, initially the number of nearest neighbors is rapidly reduced as the clusters form plates and after z 100 » the reduction rate diminishes as the plates are coarsening. For the γ phase the process is quite different. The amount of nearest neighbors reduces gradually as the system starts to grow in the z direction and the columns bend between the alpha and beta plates. Then, the average number of nearest neighbors for the phase decreases also slowly as the columns coarse.

Condensation of the mixture parameters
From the phase-field model the average thermal conductivity of the solid phase, the diffusion coefficient in the liquid phase and the average spacing between the phases can be extracted.
Since the solid structure formed during the solidification process is different in the growth direction from that formed in the directions perpendicular to the growth one, the solid conductivity is anisotropic. For carrying out the homogenization for the calculated solid field, a fixed temperature difference is imposed in the normal planes to a given direction, while adiabatic boundary conditions are imposed in the other two directions. Figure 15 shows an example of the configuration for the homogenization procedure in the z-direction, while Dirichlet boundary conditions are imposed at z 0 * = and z L z * * = by setting a fixed temperature difference between these two planes, on the planes normal the x and y directions adiabatic boundary conditions are used. The temperature field is solved for this problem according to equation (40). The equivalent thermal conductivity in the z-direction is then numerically computed as,  Average number of nearest neighbors for the cross section of the system as a function of z*. The number of nearest neighbors for a phase is the number of independent phases of the same kind that lies within a circle of radius R 1 * = in the cross section. The average is calculated over all the independent phases of a kind over the cross section. The oscillatory behavior in the number of nearest neighbors indicates an oscillatory behavior in the solid phases of the system in the z direction.
The thermal conductivity in the other two directions is numerically calculated in a similar form. Finally, the thermal conductivity of the system is expressed as a second rank tensor in the following form, In a similar form, the diffusion coefficient in the liquid next to the interface is homogenized for maintaining the solute flux. For doing so, the system is let to evolve until the it reaches a stationary state and then the diffusion coefficient is homogenized as, Finally, the separation between the formed solid phases is needed in order to know the separation of the columns that may penetrate into the solid-liquid field during the solidification process. The separation of each phase is calculated as the average distance between the centers of phases of the same type for the plane z L z * * = . The average separation of the system is then computed as the sum of the separation of the phases of the same type weighted with the area of the phase. The obtained value for the separation between phases is 3.5 sp L l = .

Conclusion
In the first section of the present work, a semi-analytical model for the solidification of the ternary eutectic FLiNaK salt has been developed for steady growth conditions. The distribution of the solid phases formed during the solidification process has been determined by minimizing the energetic cost during this process (interface undercooling). It was found that the solidified structure consists of thin long plates of the LiF and KF phases with the NaF phase present in the form of columns between these plates. By analyzing the undercooling of the solidification interface, it was observed that the energy implied in solute migration is approximately 2 times larger than the energy necessary for interfaces formation during the solidification process in the FLiNaK system. In the second part of the present work, the solidification of eutectic molten FLiNaK was studied in the mesoscale by using a phase-field model, solved by means of a pseudo-spectral method in a representative domain. It was observed that when the solidification process reaches steady state, the solid structure computed closely resembles the one previously obtained with the semi-analytical model. During the evolution towards steady solidification, the LiF and KF phases start growing as columns and then merge in order to form plate-like structures. The numerical simulations also show that the NaF phase starts growing in thin columns and, as solidification progresses, these columns bend towards the junctions of the plates of the other two phases. Afterwards, these columns cluster parallel to these junctions, growing as columns of a bigger mean diameter. During the growing process, it is observed that the interaction between forces involved in solute redistribution in the liquid and in the solid-solid interface generation, due to surface tension, produces oscillatory patterns in the growth of the three phases.
Finally, the thermal conductivity tensor, the solute diffusion coefficient for the liquid and the spacing between the phases were numerically computed from the calculated solid microstructure. A highly anisotropic thermal conductivity have been observed in the solid, being approximately 8 times larger in the growth direction of the solid than in the perpendicular directions. The condensed parameters can be introduced in a computationally less expensive macroscale mixture model. This will allow to compute the large scale solidification scenarios for molten salts found in industrial applications and, in particular for the present work, for the design of different components in the MSFR.
As a future perspective, fluid convection caused by temperature and solute natural convection needs to be considered. In addition, forced convection may be introduced in the solidification process to generate solidification conditions closer to the ones existing in the nuclear and solar concentration industries. The present work is therefore intended to be a first step toward the development of realistic and accurate solidification models for molten salts.