# Phase-field modeling for freeze crystallization in a binary system

Xiaoqian Huang, Aurélie Galfré, Françoise Couenne, Claudia Cogné

► To cite this version:

Xiaoqian Huang, Aurélie Galfré, Françoise Couenne, Claudia Cogné. Phase-field modeling for freeze crystallization in a binary system. Proceedings of the 34 th European Symposium on Computer Aided Process Engineering / 15 th International Symposium on Process Systems Engineering (ES-CAPE34/PSE24), Jun 2024, Florence, Italy. hal-04928978

HAL Id: hal-04928978

https://hal.science/hal-04928978v1

Submitted on 20 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Phase-field modeling for freeze crystallization in a binary system

Xiaoqian Huangª, Aurélie Galfréª, Françoise Couenneª, Claudia Cognéª

ªUniversité Claude Bernard Lyon 1, LAGEPP UMR5007 CNRS, 43 boulevard du 11 Novembre 1918, Villeurbanne 69100, France
Xiaoqian.huang@univ-lyon1.fr

## Abstract

This article investigates the freeze crystallization process from water-sodium chloride solution using the Phase-field model. This process (from 273.15 K to the eutectic point 252.05 K) generates two phases: a solid phase comprising solely ice and a liquid phase consisting of a concentrated solution. The phase-field model is employed to simulate the process and depict the specific dendritic structure of the solid. The challenge is to reproduce the exclusion of the solute in the ice while respecting the thermodynamic consistency of the system. The article aims to provide a straightforward thermodynamic approach to address this complexity and reproduce the solute trapping mechanisms. Additionally, the model is extended to a more concentrated solution. The conclusion suggests that lower growth velocity is favorable for ice purification and emphasizes the importance of avoiding side dendrites.

**Keywords:** Phase-field model, thermodynamics, freeze crystallization, saline solution

## 1. Introduction

Freeze crystallization, a commonly employed separation process in the food, pharmaceutical, and water treatment industries, is known for its efficiency in preserving delicate molecules and its energy efficiency. The conventional approach for modeling this process in chemical engineering is as follows: liquid and solid phases are in two distinct regions separated by a sharp interface. Conservation equations are solved in each phase and explicit algebraic conditions are imposed at the interface (Najim and Krishnan, 2022). However, this approach encounters challenges in accurately tracking the moving interface, often relying on assumptions of a sharp and discontinuous interface. Phase-field (PF) model addresses these issues by monitoring the ice front and its complex dendritic morphology by providing a spatially diffuse interface along a finite width. This diffuse interface is represented through an order variable $\varphi$. This model serves as a widely employed numerical tool in material sciences, particularly for the investigation of microstructure evolution. While the majority have been developed to address metal solidification, a few have reported studies on ice crystallization kinetics. One of the major challenges in ice crystallization modeling is the vacancy of the solute in the solid phase. To avoid salt trapping in ice, the system's thermodynamics has to be consistent and well-adjusted in the model equations. Yuan et al. (2020) studied the ice crystallization from salt solution without addressing this issue: solute is incorporated in the solid phase, which is not thermodynamically consistent. Van der Sman (2016) and Li and Fan (2020) tried to tackle this problem by introducing an interaction coefficient of components depending on $\varphi$ to define each phase's free energy. In this article, this challenge is addressed by a straightforward thermodynamic approach. In the first part,
the PF model is introduced. Then, we focus on the introduction of a pseudo-component in the model equations to mimic thermodynamically the salt behaviour. Subsequently, the developed model will be applied to the freeze crystallization of a binary system composed of H₂O-NaCl (A-B).

## 2. Phase-field model

The system under consideration consists of a small volume of H₂O-NaCl (A-B) liquid mixture. It undergoes a phase transition between solid and liquid: pure ice crystal grows due to the undercooling degree whereas the solute (salt) remains in the liquid phase. The PF model deals with the interface-tracking problem by using a continuous order variable φ. The system comprises a diffuse interface, a solid phase and liquid phase. φ takes value 0 for the solid phase (S), 1 for the liquid phase (L), and intermediate values [0, 1] in the diffuse interface. The PF model is based on the Ginzburg-Landau free energy functional of the biphasic system, energy, mass and the so-called Allen-Cahn partial differential equations defined over the whole system.

### 2.1. Governing equations

The Landau-Ginzburg Gibbs free energy functional G (J) guarantees the thermodynamic consistency of the heterogeneous system. It reads:

$$G = \int \left[ g(T, \varphi, x_B) + \frac{\epsilon_{\varphi}^2}{2} |\nabla \varphi|^2 \right] dV \quad (1)$$

The 1st term represents the free energy density g (J/m³) (Eq. (2)) depending on φ, molar fraction of the solute x_B (-), and temperature T (K); the 2nd term is the energy gradient in the interface with the interface gradient coefficient ε_φ ((J/m)¹/²):

$$g(T, \varphi, x_B) = g^S(T, x_B^S) + h(\varphi)(g^L(T, x_B^L) - g^S(T, x_B^S)) + wp(\varphi) \quad (2)$$

g comprises each bulk energy (g^S and g^L) and the energy barrier of the phase transition w(J/m³); h(φ) = φ³(6φ² - 15φ + 10) is a monotonously increasing polynomial and p(φ) = φ²(1 - φ)² is a double-well function: they are chosen to ensure local minima with respect to each phase and guarantee the thermodynamic conditions (Callen, 1985). The molar fraction at the interface x_B is related to liquid phase x_B^L and solid phase x_B^S as:

$$x_B = x_B^S + h(\varphi)(x_B^L - x_B^S) \quad (3)$$

The bulk energy of liquid phase is obtained by applying a mixture law:

$$g^L(T, x_B^L) = x_A^L \mu_A^L + x_B^L \mu_B^L \quad (4)$$

Where, for an ideal phase:

$$\mu_i^L(T, x_B^L) = \mu_i^{L,o}(T) + RT \ln x_i^L \quad (i = A \text{ or } B) \quad (5)$$

o indicates the pure component. Similar equations are used for the solid phase. By assuming local equilibrium conditions, we have:

$$\frac{\partial g^L}{\partial x_B^L} = \frac{\partial g^S}{\partial x_B^S} \Rightarrow \mu_B^L(T, x_B^L) - \mu_A^L(T, x_B^L) = \mu_B^S(T, x_B^S) - \mu_A^S(T, x_B^S) \quad (6)$$

The dynamical evolution of φ is given by the Allen-Cahn equation:

$$\frac{\partial \varphi}{\partial t} = - M_{\varphi} \frac{\delta G}{\delta \varphi} = M_{\varphi} [ \epsilon_{\varphi}^2 \nabla^2 \varphi - h'(\varphi) \Delta g - w p'(\varphi) ] \quad (7)$$
With $M_φ$ ((J.s)$^{-1}$) the interface mobility coefficient and $\frac{\delta G}{\delta φ}$ the variational derivative of G wrt φ. The notation ′ denotes derivatives wrt φ. The thermodynamic driving force $\Delta g$ is:

$$ \Delta g = g^L(T, x_B^L) - g^S(T, x_B^S) - (x_B^L - x_B^S) \frac{\partial g^L}{\partial x_B^L} \quad (8) $$

Only diffusive transport mechanism is considered in the model. The dynamic evolutions of molar fraction $x_B$ and temperature T are given by:

$$ \frac{\partial x_B}{\partial m \partial t} = - \nabla \cdot J_B = \nabla \cdot M_B(\varphi) \nabla \frac{\delta G}{\delta x_B} \quad (9) $$

$$ \frac{\partial e}{\partial t} = - \nabla \cdot J_e = \nabla \cdot \lambda(\varphi) \nabla T \quad (10) $$

$e$ (J/m$^3$) is the internal energy density, $J_B$ (mol/m$^2$/s) and $J_e$ (W/m$^2$) are respectively the densities of solute and heat fluxes. $M_B(\varphi) = \frac{D_B x_B (1-x_B)}{v_m R T}$ (mol$^2$/J.s.m)) is the mobility coefficient of B, with $v_m$ (m$^3$/mol) the constant molar volume, $D_B(\varphi)$ (m$^2$/s) the diffusivity coefficient of B in A and $\lambda(\varphi)$ (W/(m.K)) the thermal conductivity coefficient depending on the $\varphi$ evolution.

The internal energy is postulated by:

$$ e(T, \varphi, x_B) = x_A e_A(T, \varphi) + x_B e_B(T, \varphi) \quad (11) $$

Where

$$ e_i(T, \varphi) = e_i^S(T) + h(\varphi)(e_i^L(T) - e_i^S(T)) \quad (i = A \text{ or } B) \quad (12) $$

## 2.2. Interface equilibrium conditions

The parameters such as $\epsilon_φ$, $w$ and interface thickness $\xi$ (m) are related to physical properties as interface energy $\sigma$ (J/m$^2$) by considering the equilibrium conditions at interface: $\Delta g = 0$, $\frac{\partial \varphi}{\partial t} = 0$ and at melting temperature of the solution $T^m$ (K). The interface thickness and surface tension are determined (Kim et al., 1999): $\xi = \frac{4\sqrt{2}}{\sqrt{w}} \epsilon_φ$ and $\sigma = \frac{\sqrt{2w}\epsilon_φ}{6}$.

## 2.3. Pseudo-component

The PF method requires the definition of the Gibbs energy density of each phase (Eq. (4)). At equilibrium, the thermodynamic driving force $\Delta g$ must be zero (Eq. (8)). Furthermore, the local equilibrium condition (Eq. (6)) needs to be respected. In the real phase diagram of the system (solid curves in Fig. 1 represent real phase diagram with $T^E$ and $x_B^E$, the temperature and the concentration at eutectic point for H$_2$O-NaCl mixture), the solidus curve is along the y-axis, i.e. variable $x_B^S$ values zero within the study conditions. To mimic this behavior, the pseudo-component is used and enables to satisfy the PF method requirements (Bayle, 2020).

Only two properties are needed for the determination of the pseudo-component: pseudo melting point $T_B^m$ at $x_B = 1$, and pseudo melting latent heat $L_B^m$ (J/mol). Assuming the solution is ideal and neglecting the impact of thermal capacities, the liquidus equation is obtained by equal chemical potential of each component:
$$ \mu_i^L(T) - \mu_i^S(T) = L_i^m \left(1 - \frac{T}{T_i^m}\right) + RT \ln \left(\frac{x_i^L}{x_i^S}\right) = 0 \quad (i = A \text{ or } B) \quad (13) $$

By fixing $T_B^m$ smaller than $T_A^m$ but higher than absolute zero, a pseudo latent heat $L_B^m$ constant is then estimated by minimizing the sum of the squared errors between the temperature calculated for the ideal liquidus (eq.13) with real solidus (dashed curve ① in Fig.1) and simulated results. It could be noticed that the partition coefficient $k = \frac{x_B^S}{x_B^L}$ is zero for real solidus and the real liquidus requires the activity coefficient of the solute. Only the ideal solution is considered in this paper because the system is dilute. The pseudo-component is evaluated for $x_B = [0, 0.06]$ and $T^m = [260.15 K, 273.15 K]$. The difference between the estimated liquidus and ideal liquidus is small, i.e. 0.002 K.

<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of pseudo-component properties determination</figcaption></figure>

In the encircled part of Fig. 1, the pseudo solidus is represented by dashed curve ② with the pseudo properties estimated. With the known $T_B^m$, $L_B^m$ and liquid phase properties, the other thermodynamic properties required in the PF model, such as internal energy and heat capacity, can be easily calculated.

## 2.4. Numerical implementation

<table><caption>Table 1: Physical and numerical parameters</caption><thead><tr><th>Physical properties</th><th>Phase</th><th>Values</th><th>Numerical parameters</th><th>Values</th></tr></thead><tbody><tr><td rowspan="2">Thermal conductivity coefficient λ (W/(m.K))</td><td>L</td><td>0.562</td><td>Mesh number</td><td>300</td></tr><tr><td>S</td><td>2.2190</td><td>Grid size dx=dy (µm)</td><td>0.29</td></tr><tr><td rowspan="2">Diffusion coefficient of the solute D<sub>B</sub> (m<sup>2</sup>/s)</td><td>L</td><td>6.10<sup>-9</sup></td><td>Time step dt (ns)</td><td>5.23</td></tr><tr><td>S</td><td>1.73310<sup>-13</sup></td><td>Interface thickness ξ (m)</td><td>3.2510<sup>-6</sup></td></tr><tr><td>Surface tension σ (J/m<sup>2</sup>)</td><td></td><td>0.0773</td><td>Interfacial gradient coefficient ε<sub>φ</sub> ((J/m)<sup>1/2</sup>)</td><td>3.7310<sup>-4</sup></td></tr><tr><td>Molar volume v<sub>m</sub> (m<sup>3</sup>/mol)</td><td></td><td>18.10<sup>-6</sup></td><td>Magnitude of energy barrier w (J/m<sup>3</sup>)</td><td>3.3710<sup>-7</sup></td></tr><tr><td>Melting temperature of H<sub>2</sub>O T<sub>A</sub><sup>m</sup> (K)</td><td></td><td>273.15</td><td>Strength of anisotropy δ (-)</td><td>0.08</td></tr><tr><td>Latent heat of H<sub>2</sub>O L<sub>A</sub><sup>m</sup> (J/mol)</td><td></td><td>6030</td><td>Number of branches m (-)</td><td>6</td></tr><tr><td>Pseudo melting temperature of NaCl T<sub>B</sub><sup>m</sup> (K)</td><td></td><td>2</td><td>Interface mobility coefficient M<sub>φ</sub> ((J.s)<sup>-1</sup>)</td><td>2, 4</td></tr><tr><td>Pseudo latent heat of NaCl L<sub>B</sub><sup>m</sup> (J/mol)</td><td></td><td>77.3</td><td></td><td></td></tr><tr><td>Initial liquid temperatures T<sup>0</sup>(K)</td><td></td><td>258.15, 263.15, 268.15</td><td></td><td></td></tr><tr><td>Initial liquid concentration x<sub>B</sub><sup>θ</sup> (-)</td><td></td><td>10<sup>-5</sup>, 10<sup>-4</sup></td><td></td><td></td></tr></tbody></table>

The model (Eq. (7) (9) (10)) is implemented in MATLAB with the parameters in table 1. The time integration is solved by the explicit Euler method. Discretization and divergence are respectively calculated by the finite difference method and 9-points
Laplacian. The boundary conditions are fixed to be zero Neumann conditions. The initial conditions are:

liquid phase $x^2 + y^2 \ge r^2 : \varphi = 1, T = T^0, x_B^L = x_B^0$;

solid phase $x^2 + y^2 \le r^2 : \varphi = 0, T = T_A^m, x_B^S = 0$.

with $r$ the initial radius of the nuclei, and solid phase at $T_A^m = 273.15$ K the ice melting temperature, several initial liquid temperatures $T^0$ and concentrations $x_B^0$ are studied.

The anisotropy function is included in gradient coefficient $\epsilon_\varphi$:

$$ \epsilon_\varphi(\theta) = \epsilon_\varphi^0 (1 + \delta \cos(m\theta)) \qquad (16) $$

with $\delta$ the strength of the anisotropy, $m$ the number of branches and $\theta$ the angle between the direction of normal velocity and x-axis: $\theta = \arctan(\frac{\nabla\varphi_y}{\nabla\varphi_x})$ (Kobayashi, 1993).

### 3. Results and discussion

During freeze crystallization, the solute is not incorporated in the ice and accumulates at the liquid-solid interface. If the solute's diffusion velocity is slower than the interface mobility, it becomes trapped in the solid phase, resulting in a non-zero partition coefficient by forming solute brines or pockets for experimentation. Thermodynamically, this manifests as a chemical potential jump at the interface, indicating the local solute trapping effect on crystal growth and solute distribution. Regarding the temperature field, the simulated system experiences heating due to the released latent heat, leading to a reduction in undercooling and a subsequent slowdown in growth velocity. PF model effectively reproduces these phenomena.

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Concentration profile for $\Delta T = 15K$ at $t = 1.52ms$</figcaption></figure>

<caption>Figure 2: Dendrites growth for different undercooling with $x_B^0 = 10^{-5}$ at $t = 1.52ms$</caption>

In our study, the impacts of undercooling degree ($\Delta T = T^0 - T_A^m$) and solution concentration are investigated. Higher undercooling re-
sults in increased crystal growth velocity and more robust branches (Fig. 2). The solute distribution is depicted in cross-sections along the X, Y axis, and the diagonal section (Fig. 3). Solute accumulates between the main branches with lower growth velocity, while less amount of solute is observed on the branches' tips with faster velocity. The diagonal section shows the solute integration in side branches with the second minor jump. The X-axis cross-section reveals that a lower undercooling leads to an earlier side

<figure><img src="image_5.png" /><figcaption>Figure 4: Concentration profile along the X axis for (a) different undercooling (5,10,15K) with $x_B^o = 10^{-5}$ and (b) more concentrated solution $x_B^o = 10^{-4}$ at $t = 1.52\,\text{ms}$</figcaption></figure>

branches growth due to the limit of solute accumulation (Fig.4 (a)). The solute trapping effect intensifies with increasing concentration. Figure 4 (b) demonstrates that the appropriate undercooling and mobility coefficient $M_\phi$ can alleviated spurious effects by adjusting the relative velocity of crystal growth and solute diffusion. All these results support also the applicability of our pseudo-component approach to higher concentration while maintaining ice purity.

## 4. Conclusions

This article investigated the freeze crystallization process of saline solution by PF model, introducing a novel approach with a pseudo-component to address the solute vacancy for thermodynamic equilibrium and extending the model to more concentrated solution. The simulation results successfully reproduced the solute trapping mechanisms, concluding that lower crystal growth velocity leads to better ice purification. Therefore, avoiding side branches is crucial. The Further work will focus on introducing the activity of the solution and establishing the connections between the model and experimental results.

## References

* Bayle, R., 2020. Simulation des mécanismes de changement de phase dans des mémoires PCM avec la méthode multi-champ de phase.

* Callen, H.B., 1985. Thermodynamics and an Introduction to Thermostatistics. Wiley.

* Kim, S.G., Kim, W.T., Suzuki, T., 1999. Phase-field model for binary alloys. PHYSICAL REVIEW E.

* Kobayashi, R., 1993. Modeling and numerical simulations of dendritic crystal growth. Physica D: Nonlinear Phenomena 63, 410-423. https://doi.org/10.1016/0167-2789(93)90120-P

* Li, J.-Q., Fan, T.-H., 2020. Phase-field modeling of macroscopic freezing dynamics in a cylindrical vessel. International Journal of Heat and Mass Transfer 156, 119915. https://doi.org/10.1016/j.ijheatmasstransfer.2020.119915
Phase-field modeling for freeze crystallization in a binary system

Najim, A., Krishnan, S., 2022. A similarity solution for heat transfer analysis during progressive freeze-concentration based desalination. International Journal of Thermal Sciences 172, 107328. https://doi.org/10.1016/j.ijthermalsci.2021.107328

van der Sman, R.G.M., 2016. Phase field simulations of ice crystal growth in sugar solutions. International Journal of Heat and Mass Transfer 95, 153–161. https://doi.org/10.1016/j.ijheatmasstransfer.2015.11.089

Yuan, H., Sun, K., Wang, K., Zhang, J., Zhang, Z., Zhang, L., Li, S., Li, Y., 2020. Ice crystal growth in the freezing desalination process of binary water-NaCl system. Desalination 496, 114737. https://doi.org/10.1016/j.desal.2020.114737