Integrating Calphad and finite volume method for predicting non-equilibrium solidification of lithium metasilicate

Abstract

Efficient recycling of lithium metasilicate (Li₂SiO₃) from lithium-containing slag via a pyrometallurgical route demands a comprehensive understanding of its solidification process in the slag reactor. A simulation framework is developed to predict the heterogeneous phase distribution of Li₂SiO₃, the temperature and velocity fields considering density changes in the solidifying melt, on the apparatus scale. This framework integrates thermodynamic models via calculation of phase diagrams with the enthalpy-porosity technique and the volume of fluid method within a finite volume approach, ensuring thermodynamic consistency and adherence to mass balance. Thus, the formation of Li₂SiO₃ from the liquid slag composed of Li₂O-SiO₂ is described in space and temporal fields. Thereby, the interrelationship between the temperature field, enthalpy field, velocity field, and phase distribution of Li₂SiO₃ is revealed. It is shown that the lower temperature on reactor boundaries prompts the earlier formation of Li₂SiO₃ in the vicinity of the boundaries, which subsequently induces a downward flow due to the higher density of Li₂SiO₃. The predicted global mass fraction of Li₂SiO₃ under non-equilibrium conditions is 11.5 wt % lower than that calculated using the global equilibrium assumption. This demonstrates the global non-equilibrium behavior on the process scale and its consequences on slag solidification.

Graphical abstract

Keywords

solidification / thermodynamic modeling / volume of fluid / finite volume method / process simulation / Li ₂SiO ₃ crystallization

1 Introduction

Nowadays, lithium metasilicate (Li₂SiO₃), as a key component in the pseudo-binary system Li₂O-SiO₂, is attracting increasing attention in the application of electrodes in lithium-ion batteries (LIBs) [1−10]. Due to its fast Li⁺ ion conductivity property with a three-dimensional Li⁺ path, Li₂SiO₃ can be mixed with LiNi_0.5Mn_1.5O₄ to form new cathode composites so as to improve cycling stability, rate performance, and lower charge transfer resistance [1]. Recently, novel double-buffer-phase embedded nanocomposite Si/TiSi₂/Li₂SiO₃ has been proposed, and it shows good thermal stability and greater rate capability at evaluated temperatures [3]. Furthermore, Li₂SiO₃ is also widely applied as a coating material to greatly enhance the Li⁺ diffusion rate [7,9,11], the electrode structural stability [8], and electrochemical performance [12]. Other uses of Li₂SiO₃ also include serving as silicate-based sorbents for CO₂ capture [13,14], acting as a heterogeneous catalyst for producing biodiesel from soybean oil [15], reducing the sintering temperature in the low-temperature co-fired ceramic technology [16], and serving as tritium breeder material in fusion reactors [17]. The wide usage of Li₂SiO₃ also adds to the increasing demand for lithium, which is categorized as a critical raw material [18]. Therefore, ensuring lithium availability in the supply chain is crucial. However, the current processing rate may not fulfil the market demand in terms of lithium and its compounds [19], highlighting the need for energy-efficient recycling processes from industrial waste streams, e.g., via tailored slags from pyrometallurgical recycling routes of spent LIBs [20,21]. In our previous works [22−24], equilibrium thermodynamic modeling of slag solidification was performed under the framework of calculation of phase diagrams (Calphad) [25−29]. This method is a powerful tool for the development of materials and processes, facilitating the Calphad and the corresponding thermodynamic properties, such as enthalpy and latent heat, for high-entropy alloys and oxide systems [22,24,30,31]. In recent years, this method has been extensively applied to the targeted design of slag phase compositions, enabling the enrichment and optimization of vital elements within a single target phase, such as enrichment of lithium in the target phase γ-LiAlO₂ in the context of engineering artificial minerals of LIBs recycling. This is important for making downstream processes such as mechanical separation and froth flotation more efficient. Since solidification rate of γ-LiAlO₂ is fast and kinetic effects only occur at high cooling rates [32], thermodynamic equilibrium models via Calphad [25−29] can be successfully applied in this context, if external kinetic requirements in terms of heat extraction are fulfilled [24,33]. However, for Li₂SiO₃, it could be emphasized that internal kinetic effects during solidification occur already at lower cooling rates compared to γ-LiAlO₂. As an example, in Ref. [34,35], it was shown that for local thermal rates of 5 K·min^–1, significant internal kinetic effects should be taken into account. Moreover, on the reactor scale, in addition to the local kinetic effects, large spatial and temporal temperature fields in particular, as well as density and flow fields, play a central role in solidification. Computational fluid dynamics (CFD) [36] in combination with the volume of fluid (VOF) [37,38] and the enthalpy-porosity technique [39,40] offers approaches to take these phenomena into account. While there are numerous works using the enthalpy-porosity technique [39,40] for phase change materials (PCMs) [41−44] solidification from the energy storage sector, this methodology has been used less frequently for the solidification of slags. Moreover, in the traditional modeling framework, such as enthalpy-porosity technique [39,40] for PCMs solidification [41−44], the density is typically considered either constant or as a function of temperature. Conversely, when density is treated as temperature-dependent, the fixed-grid method [45] is generally employed. In Ref. [45], the volume of the entire fluid domain is presumed to remain constant, which precludes the ability to adhere to the overall mass balance [45−47]. In this work, stronger density changes during solidification are expected. Thus it is the goal to develop a simulation framework for slag solidification, where spatial and temporal flow, density, and temperature fields are considered globally and locally phase transition and phase evolution are predicted based on thermodynamic models, under the assumption of local thermodynamic equilibrium. The density of the mixture is considered to be temperature-dependent, and a dynamic mesh approach is employed, to ensure the maintenance of overall mass balance. The framework is applied to simulate the solidification of Li₂SiO₃ from Li₂O-SiO₂ melt. The grid size of the model is to be designed in such a way that the hypothesis of local thermodynamic equilibrium applies at the local cell level, where limits for this particular material system are calculated based on the thermodynamic extremal principle (TEP) according to Ref. [34]. Accordingly, the enthalpy change in the transformation zone is calculated also with thermodynamic models based on Calphad [25−29]. With this framework in place, the time evolution of the heterogeneity in solid phase distribution at the processing scale together with the interdependent spatial inhomogeneity of temperature fields and flow patterns is simulated. Such an investigation contributes to the understanding and predictive capability of oxide solidification in processing environments, providing crucial insights for efficient apparatus and process design and operation.

2 Methods

2.1 Model and VOF

In this study, a two-dimensional numerical model of the solidification process is developed, representing a cylindrical reaction chamber containing a liquid slag characterized by a chemical composition of 55 mol % SiO₂ and 45 mol % Li₂O. For the sake of demonstration, the diameter is set at 55 mm, and the height is at 130 mm. In this case, the reaction chamber is completely filled with liquid slag. The temperature of the entire fluid domain and boundaries was assumed to be uniform at 1455 K in the initial state. Following this, a cooling rate of 50 K·h^–1 is set at the four boundaries of the reaction chamber. To consider a multiphase system, the VOF method, according to Ref. [37,38] is considered. Specifically, as suggested by Richter et al. [47], we applied the VOF method [37,38] coupled with Calphad to describe the volume fraction of Li₂SiO₃ and liquid slag, denoted as α_s and α_l, respectively. To accurately describe the solidification process in terms of phase change based latent heat and its impact on the temperature, flow, and density field of this slag system, a reformulated enthalpy-porosity technique that derives its enthalpy information directly from the Calphad framework [25−29] is developed.

2.2 Enthalpy-porosity technique coupled with Calphad, and Calphad-based phase fraction calculation at the local scale (finite volume element)

In the enthalpy-porosity technique [39,40], enthalpy of the mixture is divided into two terms, i.e., sensible heat and latent heat content. Furthermore, the latent heat is assumed to vary linearly with temperature in the transition zone. However, in this work, the enthalpy of the mixture is directly deduced from the Gibbs energies from Calphad framework [25−29]. Thus for each individual cell of the simulation domain, at a particular composition and experiencing a certain temperature, the specific enthalpy of the mixture is directly calculated based on the Calphad database [22,48] from the Gibbs energy models of the respective phases. For the sake of demonstration, Fig.1(a) represents the variation of the specific enthalpy of the binary Li₂O-SiO₂ system with an initial composition of 55 mol % SiO₂ and 45 mol % Li₂O over temperature for a single finite volume. This shows for the particular material system a strong deviation from linearized latent heat calculations, as is traditionally considered in the enthalpy-porosity method [39,40]. The associated liquid fraction for the same system and over the same temperature range as is derived from the Calphad database [22,48] is presented in Fig.1(b). A significant deviation from the liquid fraction value calculated by the linear change assumption with temperature [39,40] is observed for this material system. The deviation is most pronounced as the system approaches the solidus temperature. The overall nonlinearity of the liquid fraction might be attributed to two factors: first, with the initial composition of 55 mol % SiO₂ and 45 mol % Li₂O, the phase boundary at liquidus temperature is parabolic. Second, the eutectic reaction at solidus temperature is present:

(1)\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}\leftrightarrow {\mathrm{L}\mathrm{i}}_2{\mathrm{S}\mathrm{i}\mathrm{O}}_3+{\mathrm{L}\mathrm{i}}_2{\mathrm{S}\mathrm{i}}_2{\mathrm{O}}_5-\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}.

When the system reaches the solidus temperature (Fig.1(a)), the liquid phase undergoes a sudden transformation into Li₂Si₂O₅-high temperature form and Li₂SiO₃. In this work, Calphad framework [25−29] is applied to each local finite volume cell to calculate the specific enthalpy and phase fraction via polynomial mapping at the local temperature. Then, the calculated local volume fraction of phases denoted by α_l, and α_s is provided as an input information to VOF [37,38] and the local specific enthalpy is linked to the reformulated enthalpy-porosity technique, such that the flow, temperature, and density fields can be globally simulated. Besides, for the liquid slag, liquid fraction β_l is set to 1, and β_s to 0 for the solid phase. A schematic illustration of these concepts and the geometry of the reaction chamber, which is applied for the process simulation is presented in the Fig.2. The upper boundary of the reaction chamber is configured to be movable in order to maintain mass balance as the density varies during the solidification process.

2.3 Global scale and connecting to local information

The VOF [37,38] based CFD that simulates the temperature and flow fields in this work for a process scale slag solidification, utilizes the cell level solidification information in terms of phase fraction and enthalpy, that are calculated by the coupled framework of Calphad [25−29] and enthalpy-porosity technique [39,40]. Herein, the solid deformation of the material is not considered. Additionally, due to the extremely slow-moving velocity appearing in the system, it is assumed that no turbulence is present. Under these conditions, the mass conservation equation, energy equation and momentum equation are formulated respectively as follows [40,50]:

(2)\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left(\rho \overrightarrow{v}\right)=0,

(3)\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\rho H\right)+\mathrm{\nabla }\cdot \left(\rho \overrightarrow{v}H\right)=\mathrm{\nabla }\cdot \left(k\mathrm{\nabla }T\right)+S_H,

(4)\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left(\rho \overrightarrow{v}\right)+\mathrm{\nabla }\cdot \left(\rho \overrightarrow{v}\overrightarrow{v}\right)=-\mathrm{\nabla }p+\mathrm{\nabla }\cdot \left(\overline{\overline{\tau }}\right)+\rho \overrightarrow{g}+F.

In the above equations, p is pressure, \rho represents density, t represents time, \overrightarrow{v} is velocity vector, H is enthalpy, k is heat conductivity, S_H is the source term, which is set to zero for the particular case studied in this work. \overrightarrow{g} is gravitational acceleration, \overline{\overline{\tau }} is stress tensor, which can be expressed as according to [50]:

(5)\overline{\overline{\tau }}=\mu \left[\left(\mathrm{\nabla }\overrightarrow{v}+\mathrm{\nabla }{\overrightarrow{v}}^T\right)-\frac{2}{3}\mathrm{\nabla }\cdot \overrightarrow{v}I\right].

Here, {\overrightarrow{v}}^T is the transpose of the velocity vector, I is the unit tensor. In Eq. (5), besides the gravity force, there is another body force term F for modifying the momentum in the mush region, which can be described according to [40,47,51]:

(6)F=\frac{{\left(1-\beta \right)}^2}{{\beta }^3+\epsilon }{\mathit{A}}_{\mathrm{m}\mathrm{u}\mathrm{s}\mathrm{h}}\overrightarrow{v}.

In Eq. (6), \epsilon is a constant, which is introduced to avoid division by zero. A_mush is an adjustable mush zone constant. It is worth emphasizing that the liquid fraction of the mixture, β, in the present work does not follow the linear temperature relationship from the classical enthalpy-porosity technique [39,40], and is derived based on Calphad method [25−29]. The parameter β can be expressed as:

(7)\beta ={\alpha }_{\mathrm{s}}{\beta }_{\mathrm{s}}+{\alpha }_{\mathrm{l}}{\beta }_{\mathrm{l}}.

Under the condition that no material is being pulled out of the domain, the introduction of the momentum source term accounts for the strong decrease of the flow velocity of the mixture with decreasing liquid fraction. Especially when the liquid fraction is approaching zero, the predicted velocity will be close to zero.

2.4 Numerical framework

The model is implemented in Ansys fluent [52]. To discretize the geometry of the exemplified reaction chamber, different structured and quadrilateral mesh sizes, such as 0.4, 0.3, 0.2, and 0.1 mm, are applied to assess the mesh sensitivity of the numerical simulations. The details will be discussed in the results section. Regarding the solution methods employed in this study, the SIMPLE algorithm [53] was utilized for pressure-velocity coupling. For spatial discretization, the least squares cell-based scheme was adopted for gradient evaluation, while the PRESTO! Scheme [54] was applied for pressure interpolation. The momentum and energy equations were solved using the second-order upwind scheme [55]. To determine the volume fraction, the Geo-Reconstruct scheme [56] was utilized. Additionally, a first-order implicit formulation was applied for transient calculations, with a constant time step set at 1 s. Additionally, the upper boundary is modeled as a moving boundary, with its position dynamically adjusted according to a governing equation designed to maintain the overall mass balance. For this case, the dynamic mesh technique [57] was employed where layering was specifically applied to the upper boundary. Simultaneously, the left boundary, right boundary, and the interior fluid were configured to deform within the dynamic mesh settings. Surface forces and further body forces are neglected in this work. Furthermore, the applied material properties and simulation conditions in terms of initial temperature and cooling rate applied to the system are summarized in Tab.1.

In Tab.1, the volume fraction and specific heat capacity of each phase are temperature-dependent. They are calculated based on the Calphad framework [25−29] under local equilibrium conditions and mapped via polynomial functions, such that the information is transferred via a user-defined function in Ansys fluent [52].

2.5 Identifying maximal local cooling rates to ensure validation of the local thermodynamic equilibrium hypothesis based on the TEP

The validity of the local equilibrium assumption for melt solidification that has been considered in this work, depends on the external cooling rates, the transport kinetics of the material system under investigation, as well as the reaction kinetics of the solidifying phases. Hence, it is essential to elucidate the dependence of different internal kinetic effects during solidification, on the external processing conditions.

A thermodynamically consistent model based on the TEP has been developed in our previous works [34,35,64]. Within the model framework, the internal kinetic effects, as well as the migration of the solid/liquid interface, have been considered to dissipate the system Gibbs energy. Hence, first, the rate of this Gibbs energy \dot{G}, was formulated in terms of certain kinetic variables, such as the mass flux and the interface velocity that are representatives of the internal kinetic processes on the one hand and characterizes the system on the other. Next, based on the linear relation between fluxes and forces for each of the kinetic processes, the dissipation function Q was represented as a quadratic function of the same kinetic variables. Finally, considering the constraint of \dot{G}=-Q stemming from irreversible thermodynamics, a maximization of Q with respect to the kinetic variables yielded the evolution equations, hence the time evolution of the internal material and thermal profiles for a solidifying melt in dependence of the external thermal profiles could be calculated. For a detailed derivation of the evolution equations, the readers might refer to Ref. [65,66].

The TEP framework, when implemented for the Li₂O-SiO₂ system, could describe the non-equilibrium nature of Li₂SiO₃ solidification for the fast cooling of the binary melt in Ref. [34]. The model parameters could be successfully determined to capture the experimentally observed deviation of the solidified phase fraction from its equilibrium behavior.

In this work, the TEP model, with its experimentally determined parameters has been implemented inversely, such that a limiting external cooling rate can be determined below which, minimal deviation from equilibrium is predicted and the local equilibrium assumption is substantiated and chosen cell size is approved for the particularly simulated cases in this work.

3 Results and discussion

The simulation framework developed in this work has been implemented to study the solidification of Li₂SiO₃ from a liquid slag with a bulk composition of 55 mol % SiO₂ and 45 mol % Li₂O. Hence, the initial temperature of the fluid and boundaries is set at 1455 K, which is the liquidus temperature for this composition. A boundary cooling rate of 50 K·h^–1 is applied and the solidification is simulated for a duration of 1000 s, to put the particular focus on the early stage crystallization regime in this work. The time step is chosen to be 1 s, due to the very high viscosity of the slag system.

Furthermore, the simulation domain had to be discretized using a mesh of a certain dimension. For this, a sensitivity analysis was first performed to observe the influence of mesh size on the simulated solidified fractions. A structured and quadrilateral cell based grid is used for meshing. Based on the present numerical simulation, the overall mass fraction of the Li₂SiO₃ can be updated in every time step. With a coarse mesh size, such as 0.4 mm, the time-dependent phase evolution of Li₂SiO₃ exhibited frequent numerical fluctuations. A sharp increase of the global mass fraction of Li₂SiO₃ from 0 to 14.2 wt % followed by a steep decline to 7.3 wt % within only the first 200 s suggests a strong numerical instability. With the continuous refinement of the mesh, the temporal numerical fluctuation of the generated phase amounts gradually disappeared. Specifically, when the mesh size was reduced to 0.1 mm, no mass fraction jump was observed throughout the entire simulation. Therefore, in this work, the mesh size of 0.1 mm is selected for the numerical simulations.

Next, a comparison was made between the local cooling rates experienced by several representative locations within the reactor, when the boundary cooling was set at 50 K·h^–1, and the TEP determined cooling rate of approximately 10.2 K·h^–1, for which the internal kinetic effects are negligible. It was observed that the local temporal thermal gradients were in the range or even below the critical value. Hence, the assumption of the local thermodynamic equilibrium for the externally applied kinetic condition was validated.

Consequently for any boundary cooling rate below 50 K·h^–1, some of which are further discussed later in this section, the presently developed non-equilibrium model with the local equilibrium hypothesis is applicable to simulate time-dependent solidification behavior. For the cases, when the reactor is subjected to higher cooling rates at the boundary, the local thermal gradients need to be further evaluated to determine whether it falls within the TEP defined critical cooling rate. If not, this would signify that even locally, the kinetics of diffusion and phase generation are not fast enough compared to the process time scales to attain local equilibrium. This would trigger spatial heterogeneity of material and thermal profiles even at the local cell level and the corresponding solidification behavior can no longer be described by the local equilibrium hypothesis. In such cases, modeling strategies that capture the local kinetics need to be further integrated to the present framework. However, for this study, ensuring adherence to local equilibrium conditions, model calculations were first performed under the boundary cooling rate of 50 K·h^–1, with the simulation space being discretized using a 0.1 mm mesh, to elucidate the interrelationship between the temperature field, enthalpy field, flow field, and phase distribution of the Li₂SiO₃ as the liquid slag solidifies.

Fig.3 represents the heterogeneity of the temperature, velocity, and density fields formed at the reactor scale when the slag system has been cooled for 1000 s. As can be observed from Fig.3(a) the temperature on the four boundaries is generally lower than the rest of the region of the reaction chamber.

According to the numerical simulation, when the boundary temperature decreases to 1441 K at a cooling rate of 50 K·h^–1, the central region of the slag retains its initial temperature of 1455 K. This thermal lag effect thus results in a reduced generation of Li₂SiO₃, in comparison to the phase amount formed under the assumption of the global homogeneous temperature distribution. Moreover, even with identical cooling rates and initial temperatures on the four boundaries, after 1000 s, the temperature distribution near the four boundaries is not identical. This phenomenon can be attributed to two main factors. On the one hand, during the solidification process, the overall density of the mixture increases because the density of the solid phase Li₂SiO₃ is higher than that of the liquid slag. As the temperature decreases, an increasing amount of Li₂SiO₃ forms, resulting in a shrinkage of the entire fluid domain. This shrinkage causes the upper boundary to move downward. Consequently, the temperature gradient from the upper boundary to the center region becomes more pronounced than the gradient from the bottom boundary to the center region. On the other hand, the movement of the mixture resulting from density variations, as depicted in Fig.3(c) during the solidification process affects the temperature and velocity fields. The details of the velocity field are illustrated in Fig.3(b). As the temperature continues to decrease, an increasing amount of the liquid phase near the boundary transforms into the solid phase Li₂SiO₃. Additionally, due to gravity, the denser mixture near the two vertical boundaries tends to move downward. Furthermore, the downward movement is accelerated by the lowering of the upper boundary. Subsequently, the mixture from both vertical sides reaches the bottom boundary and flows in a reverse direction upward. Additionally, vortices can be observed in the flow pattern represented in Fig.3(b). The highest velocity is calculated to be around 1.4 × 10⁻⁵ m·s^–1 in the reverse flow region at the bottom area. In contrast, the movement in the center area remains very low, as the temperature in this region stays at 1455 K, preventing the production of Li₂SiO₃. Consequently, the density in this region remains constant, equal to the liquid slag density of 2025 kg·m^–3. In this work, the specific enthalpy is directly informed by the Calphad [25−29], which inherently includes the latent heat term. As temperature decreases, the liquid slag transforms into solid phase Li₂SiO₃. Due to the release of the latent heat, the temperature at local points can increase during the solidification process and lead to the remelting phenomenon. This local temperature increase also affects the phase distribution of Li₂SiO₃, as its local phase fraction, is calculated using Calphad method [25−29]. The spatial variation in the solid phase distribution further influences the density field and impacts the flow pattern. Therefore, in this numerical simulation, accurately describing the specific enthalpy of mixture is crucial for the integrated enthalpy-porosity technique [39,40]. Accordingly, the calculated specific enthalpy of the simulation domain after the system has been cooled for 1000 s is presented in Fig.4.

The absolute and specific enthalpy of the mixture calculated using the Calphad model [25−29] within the temperature range of 1455–1441 K, varies from approximately –1.550 × 10⁷ to –1.561 × 10⁷ J·kg^–1. This aligns well with the calculated enthalpy field of the mixture shown in Fig.4. Thus, this confirms that the specific enthalpy from the numerical simulation at each local position matches the specific enthalpy predicted by the Calphad method [25−29], validating the successful integration of Calphad [25−29] and enthalpy-porosity technique [39,40]. The evolution of the phase distribution of Li₂SiO₃ with time has been represented in Fig.5.

In the early stage of the solidification process, as depicted in Fig.5(a), Li₂SiO₃ is initially generated at the boundaries. However, compared to the bottom boundary, Li₂SiO₃ is less concentrated near the top boundary. As the temperature continues to decrease, this phenomenon becomes increasingly pronounced and can be attributed to two main factors. First, due to the gravitational field, the generated Li₂SiO₃, being denser than the liquid slag, tends to accumulate at the bottom boundary. The velocity field shown in Fig.3(b) also indicates downward movement along the top boundary and the presence of vortices at the bottom boundary. Additionally, the downward movement of the top boundary increases the vertical temperature gradient in this area. Compared to the initial height of 130 mm, the final height after 1000 s of simulation time is reduced to 129 mm. As illustrated in Fig.3(a), this change leads to a higher average temperature near the top boundary, which subsequently inhibits the formation of Li₂SiO₃.

The effect of the external processing conditions in terms of heat extraction-based boundary cooling rates on the solidification behavior of the slag is elucidated further in Fig.6 by comparing the generated mass fraction of Li₂SiO₃ within the different simulation times of slag solidification under different cooling conditions. The cases (i), (ii), (iii), (iv), (v), and (vi) are for boundary cooling rates of 50, 45, 40, 35, 30, and 25 K·h^–1, respectively, considering heterogeneous temperature and velocity fields. (vii) is for the sake of comparison a calculated case with a boundary cooling rate of 50 K·h^–1 with a homogeneous temperature field within the reactor chamber.

For cases (i) and (vi) in Fig.6(a) it can be observed that decreasing the cooling rate delays the formation of the first crystalline solid phase Li₂SiO₃ almost twofold, from 74 to 146 s. Moreover, since the formation rate of Li₂SiO₃ is slower at 25 K·h^–1 compared to 50 K·h^–1, only 1.4 wt % of Li₂SiO₃ is produced after 1000 s, which is less than half of the amount formed at the higher cooling rate of 50 K·h^–1. In the case of a cooling rate of 40 K·h^–1, the first crystalline phase forms at 92 s, with 2.6 wt % of Li₂SiO₃ produced after 1000 s. This amount is lower than that generated at the cooling rate of 50 K·h^–1 but higher than those at the rates which lower than 40 K·h^–1. The influence of higher density of solid and lowering of the upper boundary, which can now be precisely described under mass balance conditions, leads to the early concentration of solid phase near the bottom boundary. Furthermore, compared to the generated Li₂SiO₃ mass fraction for case (vii) denoted by the brown line in Fig.6(a), which is determined to be 14.9 wt % after 1000 s, the calculated mass fractions for cases (i), (ii), (iii), (iv), (v), and (vi) after 1000 s are significantly smaller. This stems from the fact that under the assumptions of global thermodynamic equilibrium, spatial and temporal fields have not been considered for the case (vii). Whereas for the cases (i), (ii), (iii), (iv), (v), and (vi) that better represent the conditions for a practical reaction chamber or slag reactor, the temperature and chemical composition differ at each position and hence represents non-equilibrium cases at the global scale, while locally equilibrium is maintained. Fig.6(b) illustrates the effect of varying cooling rates on the phase formation of Li₂SiO₃ across different simulation times. At each simulation time, ranging from 500 to 1000 s, the effect of varying cooling rates on the phase formation of Li₂SiO₃ appears to be approximately linear. As the cooling rate increases, the predicted phase amount of Li₂SiO₃ also rises. For instance, at the simulation time of 500 s, the mass fraction of Li₂SiO₃ increases from 0.3 to 1.0 wt % when the cooling rate is elevated from 25 to 50 K·h^–1. Moreover, the slope of the lines progressively increases between 500 and 1000 s. At 1000 s of simulation time, an additional 2.0 wt % of Li₂SiO₃ forms when the cooling rate is increased from 25 to 50 K·h^–1, compared to an increase of only 0.7 wt % at 500 s. This indicates that higher cooling rates have a more significant effect on promoting the phase formation of Li₂SiO₃ during the later stages of phase evolution.

4 Conclusions

The objective of this work is to develop a numerical simulation framework that integrates thermodynamic models via Calphad with the enthalpy-porosity technique and the VOF method within a finite volume approach. The aim was to predict phase distribution and temperature profiles during solidification processes, particularly for simulating solidification of slags on processing scale. The framework applies Calphad to integrate the enthalpy of the mixture and the temperature-dependent phase behavior in each local cell, where the local equilibrium hypothesis is applied and approved via TEP. Considering the thermal lag effect, a temperature difference of 14 K is observed between the boundary and the core of the fluid in the herein chosen simulation geometry at a cooling rate of 50 K·h^–1 after 1000 s. Due to the lower temperature near the boundaries, the generated Li₂SiO₃ is primarily distributed in the boundary area, while the central area remains in the liquid state for the observed time interval. At the same time, the density of the fluid near the boundaries increases from 2025 to approximately 2095 kg·m^–3. Owing to the significant density changes during solidification, the upper boundary which accounts for maintaining mass balance moves downward. This leads to an increase in the temperature gradient and simultaneously accelerates the flow velocity of the fluid near the top boundary. Once the fluid reaches the bottom boundary, it reverses direction and flows upward, forming vortices with a maximum predicted velocity of approximately 1.4 × 10^–5 m·s^–1. With this approach, the time- and spatial-dependent phase distribution of Li₂SiO₃ and the heterogeneous temperature field, and the corresponding flow pattern during solidification is calculated. It turned out, for the particular system studied that due to non-equilibrium and flow conditions, the mass fraction of Li₂SiO₃ is 11.5 wt % below the global equilibrium prediction. This emphasizes the necessity of taking into account heterogeneous fields and non-equilibrium on the slag processing scale, such that reactor conditions can be simulated and optimized.

References

[1] Deng Y , Mou J , Wu H , Jiang N , Zheng Q , Lam K H , Xu C , Lin D . A superior Li₂SiO₃-composited LiNi_0.5Mn_1.5O₄ cathode for high-voltage and high-performance lithium-ion batteries. Electrochimica Acta, 2017, 235: 19–31

[2] Khalifa H , El-Safty S A , Reda A , Shenashen M A , Eid A I . Anisotropic alignments of hierarchical Li₂SiO₃/TiO₂@nano-C anode//LiMnPO₄@nano-C cathode architectures for full-cell lithium-ion battery. National Science Review, 2020, 7(5): 863–880

[3] Jeong W J , Chung D J , Youn D , Kim N G , Kim H . Double-buffer-phase embedded Si/TiSi₂/Li₂SiO₃ nanocomposite lithium storage materials by phase-selective reaction of SiO with metal hydrides. Energy Storage Materials, 2022, 50: 740–750

[4] Yang S , Wang Q , Miao J , Zhang J , Zhang D , Chen Y , Yang H . Synthesis of graphene supported Li₂SiO₃ as a high performance anode material for lithium-ion batteries. Applied Surface Science, 2018, 444: 522–529

[5] Hu G , Zhang M , Wu L , Peng Z , Du K , Cao Y . Effects of Li₂SiO₃ coating on the performance of LiNi_0.5Co_0.2Mn_0.3O₂ cathode material for lithium ion batteries. Journal of Alloys and Compounds, 2017, 690: 589–597

[6] Li Y , Zhang D , Yan Y , Wang Y , Li Z , Tan X , Zhang M . Enhanced electrochemical properties of SiO₂-Li₂SiO₃-coated NCM811 cathodes by reducing surface residual lithium. Journal of Alloys and Compounds, 2022, 923: 166317

[7] Peng Z , Yang G , Li F , Zhu Z , Liu Z . Improving the cathode properties of Ni-rich LiNi_0.6Co_0.2Mn_0.2O₂ at high voltages under 5 C by Li₂SiO₃ coating and Si⁴⁺ doping. Journal of Alloys and Compounds, 2018, 762: 827–834

[8] Qiao Y , Hao R , Shi X , Li Y , Wang Y , Zhang Y , Tang C , Li G , Wang G , Liu J . . Improving the cycling stability of LiNi_0.8Co_0.1Mn_0.1O₂ by enhancing the structural integrity via synchronous Li₂SiO₃ coating. ACS Applied Energy Materials, 2022, 5(4): 4885–4892

[9] Zhao E , Liu X , Zhao H , Xiao X , Hu Z . Ion conducting Li₂SiO₃-coated lithium-rich layered oxide exhibiting high rate capability and low polarization. Chemical Communications, 2015, 51(44): 9093–9096

[10] Liu X , Ouyang B , Hao R , Liu P , Fan X , Zhang M , Pan M , Liu W , Liu K . Li₂SiO₃ modification of C/LiFe_0.5Mn_0.5PO₄ for high performance lithium-ion batteries. ChemElectroChem, 2022, 9(16): e202200609

[11] Wang D , Zhang X , Xiao R , Lu X , Li Y , Xu T , Pan D , Hu Y S , Bai Y . Electrochemical performance of Li-rich Li[Li_0.2Mn_0.56Ni_0.17Co_0.07]O₂ cathode stabilized by metastable Li₂SiO₃ surface modification for advanced Li-ion batteries. Electrochimica Acta, 2018, 265: 244–253

[12] Bai X , Li T , Dang Z , Qi Y X , Lun N , Bai Y J . Ionic conductor of Li₂SiO₃ as an effective dual-functional modifier to optimize the electrochemical performance of Li₄Ti₅O₁₂ for high-performance Li-ion batteries. ACS Applied Materials & Interfaces, 2017, 9(2): 1426–1436

[13] Kwon Y M , Chae H J , Cho M S , Park Y K , Seo H M , Lee S C , Kim J C . Effect of a Li₂SiO₃ phase in lithium silicate-based sorbents for CO₂ capture at high temperatures. Separation and Purification Technology, 2019, 214: 104–110

[14] Ortiz-Landeros J , Gómez-Yáñez C , Pfeiffer H . Surfactant-assisted hydrothermal crystallization of nanostructured lithium metasilicate (Li₂SiO₃) hollow spheres: II—textural analysis and CO₂-H₂O sorption evaluation. Journal of Solid State Chemistry, 2011, 184(8): 2257–2262

[15] Wang J X , Chen K T , Huang S T , Chen C C . Application of Li₂SiO₃ as a heterogeneous catalyst in the production of biodiesel from soybean oil. Chinese Chemical Letters, 2011, 22(11): 1363–1366

[16] Wang M , Zhang S , Yang Z , Li E , Tang B , Zhong C . Sintering behaviors and thermal properties of Li₂SiO₃-based ceramics for LTCC applications. Ceramics International, 2022, 48(19): 27312–27323

[17] Nishikawa Y , Oyaidzu M , Yoshikawa A , Munakata K , Okada M , Nishikawa M , Okuno K . Correlation between tritium release and thermal annealing of irradiation damage in neutron-irradiated Li₂SiO₃. Journal of Nuclear Materials, 2007, 367–370: 1371–1376

[18] DavidMLythS MLindnerRHarringtonG F. Future-Proofing Fuel Cells: Critical Raw Material Governance in Sustainable Energy. Cham: Palgrave Macmillan Cham, 2021, 15–33

[19] Wanger T C . The lithium future—resources, recycling, and the environment. Conservation Letters, 2011, 4(3): 202–206

[20] Elwert T , Strauss K , Schirmer T , Goldmann D . Phase composition of high lithium slags from the recycling of lithium ion batteries. World of Metallurgy-Erzmetall, 2012, 65(3): 163–171

[21] Sommerfeld M , Vonderstein C , Dertmann C , Klimko J , Oráč D , Miškufová A , Havlík T , Friedrich B . A combined pyro- and hydrometallurgical approach to recycle pyrolyzed lithium-ion battery black mass Part 1: production of lithium concentrates in an electric arc furnace. Metals, 2020, 10(8): 1069

[22] Li H , Qiu H , Schirmer T , Goldmann D , Fischlschweiger M . Tailoring lithium aluminate phases based on thermodynamics for an increased recycling efficiency of Li-ion batteries. ACS ES&T Engineering, 2022, 2(10): 1883–1895

[23] Li H , Ranneberg M , Fischlschweiger M . High-temperature phase behavior of Li₂O-MnO with a focus on the liquid-to-solid transition. Journal of the Minerals Metals & Materials Society, 2023, 75(12): 5796–5807

[24] Li H , Qiu H , Ranneberg M , Lucas H , Graupner T , Friedrich B , Yagmurlu B , Goldmann D , Bremer J , Fischlschweiger M . Enhancing lithium recycling efficiency in pyrometallurgical processing through thermodynamic-based optimization and design of spent lithium-ion battery slag compositions. ACS Sustainable Resource Management, 2024, 1(6): 1170–1184

[25] LukasHFriesS GSundmanB. Computational Thermodynamics: The Calphad Method. New York: Cambridge University Press, 2007, 1–46

[26] HillertM. Phase Equilibria, Phase Diagrams and Phase Transformations: Their Thermodynamic Basis. 2nd ed. New York: Cambridge University Press, 2007, 400–418

[27] Pelton A D . A general “geometric” thermodynamic model for multicomponent solutions. Calphad, 2001, 25(2): 319–328

[28] Pelton A D , Degterov S A , Eriksson G , Robelin C , Dessureault Y . The modified quasichemical model I: binary solutions. Metallurgical and Materials Transactions B, Process Metallurgy and Materials Processing Science, 2000, 31(4): 651–659

[29] Pelton A D , Chartrand P . The modified quasi-chemical model: Part II. Multicomponent solutions. Metallurgical and Materials Transactions A, Physical Metallurgy and Materials Science, 2001, 32(6): 1355–1360

[30] Olson G B , Liu Z K . Genomic materials design: CALculation of PHAse dynamics. Calphad, 2023, 82: 102590

[31] Wu M , Wang S , Huang H , Shu D , Sun B . Calphad aided eutectic high-entropy alloy design. Materials Letters, 2020, 262: 127175

[32] de Abreu D A , Löffler M , Kriegel M J , Fabrichnaya O . Experimental investigation and thermodynamic modeling of the Li₂O-Al₂O₃ system. Journal of Phase Equilibria and Diffusion, 2024, 45(1): 36–55

[33] Chakrabarty S , De Abreu D A , Alhafez I A , Fabrichnaya O , Merkert N , Schnickmann A , Schirmer T , Fittschen U E A , Fischlschweiger M . Kinetics of γ-LiAlO₂ formation out of Li₂O-Al₂O₃ melt: a molecular dynamics-informed non-equilibrium thermodynamic study. Solids, 2024, 5(4): 561–579

[34] Chakrabarty S , Li H , Schirmer T , Hampel S , Fittschen U E A , Fischlschweiger M . Non-equilibrium thermodynamic modelling of cooling path dependent phase evolution of Li₂SiO₃ from Li₂O-SiO₂ melt by considering mixed kinetic phenomena and time-dependent concentration fields. Scripta Materialia, 2024, 242: 115922

[35] Chakrabarty S , Li H , Fischlschweiger M . Control of interface migration in nonequilibrium crystallization of Li₂SiO₃ from Li₂O-SiO₂ melt by spatiotemporal temperature and concentration fields. ACS Omega, 2024, 9(19): 21557–21568

[36] WendtJ F. Computational Fluid Dynamics. 3rd ed. Berlin: Springer, 2009, 275–301

[37] Hirt C W , Nichols B D . Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics, 1981, 39(1): 201–225

[38] Gopala V R , Van Wachem B G M . Volume of fluid methods for immiscible-fluid and free-surface flows. Chemical Engineering Journal, 2008, 141(1-3): 204–221

[39] Voller V R , Swaminathan C R , Thomas B G . Fixed grid techniques for phase change problems: a review. International Journal for Numerical Methods in Engineering, 1990, 30(4): 875–898

[40] Voller V R , Prakash C . A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems. International Journal of Heat and Mass Transfer, 1987, 30(8): 1709–1719

[41] Muhammad M D , Badr O , Yeung H . Validation of a CFD melting and solidification model for phase change in vertical cylinders. Numerical Heat Transfer Part A, 2015, 68(5): 501–511

[42] Dallaire J , Gosselin L . Numerical modeling of solid-liquid phase change in a closed 2D cavity with density change, elastic wall and natural convection. International Journal of Heat and Mass Transfer, 2017, 114: 903–914

[43] Samar A H , Napitupulu F H , Sitorus T B . Numerical study on melting process of phase change material as thermal energy storage. IOP Conference Series. Materials Science and Engineering, 2020, 725(1): 012051

[44] MahdiM SMahoodH BHasanA FKhadomA A. Solidification enhancement of phase change material implemented in latent heat thermal energy storage. In: Proceedings of 2nd International Conference on Materials Engineering & Science. Baghdad: AIP Publishing, 2019, 020039

[45] Martínez A , Carmona M , Cortés C , Arauzo I . Experimentally based testing of the enthalpy-porosity method for the numerical simulation of phase change of paraffin-type PCMs. Journal of Energy Storage, 2023, 69: 107876

[46] Younsi Z , Naji H . A numerical investigation of melting phase change process via the enthalpy-porosity approach: application to hydrated salts. International Communications in Heat and Mass Transfer, 2017, 86: 12–24

[47] Richter O , Turnow J , Kornev N , Hassel E . Numerical simulation of casting processes: coupled mould filling and solidification using VOF and enthalpy-porosity method. Heat and Mass Transfer, 2017, 53(6): 1957–1969

[48] Konar B , Van Ende M A , Jung I H . Critical evaluation and thermodynamic optimization of the Li-O and Li₂O-SiO₂ systems. Journal of the European Ceramic Society, 2017, 37(5): 2189–2207

[49] Hesse K F . Refinement of the crystal structure of lithium polysilicate. Acta Crystallographica. Section B, Structural Crystallography and Crystal Chemistry, 1977, 33(3): 901–902

[50] MoukalledFManganiLDarwishM. The Finite Volume Method in Computational Fluid Dynamics. Cham: Springer International Publishing, 2016, 103–135

[51] Glück Nardi V , Greß T , Tonn B , Volk W . Modelling of intermetallic layers formation during solid-liquid joining of dissimilar metallic materials. Materials Science and Engineering, 2020, 861(1): 012058

[52] Ansys Fluent, Academic Research Fluent, Version 2022 R1, 2022

[53] Roy P , Anand N K , Donzis D . A parallel multigrid finite-volume solver on a collocated grid for incompressible Navier-Stokes equations. Numerical Heat Transfer Part B, 2015, 67(5): 376–409

[54] Shukla S K , Shukla P , Ghosh P . Evaluation of numerical schemes using different simulation methods for the continuous phase modeling of cyclone separators. Advanced Powder Technology, 2011, 22(2): 209–219

[55] Shyy W , Thakur S , Wright J . Second-order upwind and central difference schemes for recirculating flow computation. AIAA Journal, 1992, 30(4): 923–932

[56] TahirF. Transient analysis of air bubble rise in stagnant water column using CFD. In: Proceedings of ICTEA. Doha: IEEE, 2018

[57] Drake R , Manoranjan V S . A method of dynamic mesh adaptation. International Journal for Numerical Methods in Engineering, 1996, 39(6): 939–949

[58] Butland A T D , Maddison R J . The specific heat of graphite: an evaluation of measurements. Journal of Nuclear Materials, 1973, 49(1): 45–56

[59] Jiao J , Grorud B , Sindland C , Safarian J , Tang K , Sellevoll K , Tangstad M . The use of eutectic Fe-Si-B alloy as a phase change material in thermal energy storage systems. Materials, 2019, 12(14): 2312

[60] Sajid M , Bai C , Aamir M , You Z , Yan Z , Lv X . Understanding the structure and structural effects on the properties of blast furnace slag (BFS). ISIJ International, 2019, 59(7): 1153–1166

[61] HayashiMIshiiHSusaMFukuyamaHNagataK. Hierarchy of thermal conductivities for alkali silicates in terms of ionicity of oxygen. In: Proceedings of Internal Conference on Molten Slags, Fluxes, and Salts. Stockholm: Division of Metallurgy, KTH, Sweden, 2000

[62] Wang M , Zhong C , Yang Z , Yang H , Cao L , Qin T , Zhang S . Thermal and microwave dielectric properties of Li-Si-based ceramics. Ceramics International, 2021, 47(12): 17693–17701

[63] Bale C W , Bélisle E , Chartrand P , Decterov S A , Eriksson G , Gheribi A E , Hack K , Jung I H , Kang Y B , Melançon J . . FactSage thermochemical software and databases, 2010–2016. Calphad, 2016, 54: 35–53

[64] Chakrabarty S , Li H , Fischlschweiger M . Calphad-informed thermodynamic non-equilibrium simulation of non-isothermal solid-state reactions of magnesium aluminate spinel based on the thermodynamic extremal principle. Materialia, 2023, 28: 101723

[65] Svoboda J , Gamsjäger E , Fischer F D , Fratzl P . Application of the thermodynamic extremal principle to the diffusional phase transformations. Acta Materialia, 2004, 52(4): 959–967

[66] Abart R , Svoboda J , Jeřabek P , Povoden-Karadeniz E , Habler G . Interlayer growth kinetics of a binary solid-solution based on the thermodynamic extremal principle: application to the formation of spinel at periclase-corundum contacts. American Journal of Science, 2016, 316(4): 309–328

RIGHTS & PERMISSIONS

The Author(s) 2025. This article is published with open access at link.springer.com and journal.hep.com.cn