# An implicit regularized enthalpy Lattice Boltzmann Method for the Stefan problem

Francky Luddens, Corentin Lothodé, Ionut Danaila

► To cite this version:

Francky Luddens, Corentin Lothodé, Ionut Danaila. An implicit regularized enthalpy Lattice Boltzmann Method for the Stefan problem. 2025. hal-05070487

HAL Id: hal-05070487

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

Preprint submitted on 16 May 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
# An implicit regularized enthalpy Lattice Boltzmann Method for the Stefan problem

Francky Luddens<sup>a,*</sup>, Corentin Lothodé<sup>b</sup>, Ionut Danaila<sup>a</sup>

## Abstract

Solving the Stefan problem, also referred as the heat conduction problem with phase change, is a necessary step to solve phase change problems with convection. In this article, we are interested in using the Lattice Boltzmann Method (LBM) to solve the Stefan problem using a regularized total enthalpy model. The liquid fraction is treated as a non-linear source/sink term, that involves the time derivative of the solution. The resulting non-linear system is solved using a Newton algorithm. By conserving the locality of the problem, this method is highly scalable, while keeping a high accuracy. The newly developed scheme is analyzed theoretically through a Chapman-Enskog expansion and illustrated numerically with 1D and 2D benchmarks.

**Keywords:** Phase change, Stefan problem, Lattice Boltzmann method

**License:** This work is licensed under https://creativecommons.org/licenses/by/4.0/ (CC BY)

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

## 1. Introduction

In this article, we consider the classical two-phase Stefan problem [1, 2]. It describes the melting or solidification of a material taking into account the conduction as the main heat transfer mechanism. Numerical solvers for liquid-solid phase change systems, such as those occurring in metal casting, energy storage, cryosurgery or food freezing, are usually first benchmarked against the analytical solution of the Stefan problem [3].

The Stefan problem is well known in the mathematical community, with many theoretical results published, see for example [4] for a review. In the context of numerical simulation, several frameworks were developed to solve this problem [5]. In the finite element community, different approaches were suggested, from enthalpy based methods (e.g. [6]) to front tracking methods (e.g. [7]). For other numerical frameworks, similar approaches are developed, see for example [8] for finite volumes and [9] for finite differences.

The Lattice Boltzmann Method (LBM) was first adapted in [10] to tackle heat diffusion. Later, a chemical reaction analogy and two types of quasi particles were introduced in [11]. The resulting method involves solving a linear system of equations that was subsequently simplified in [12] using only one type of quasi particle and a phase field method. The position of the interface is recovered through the liquid fraction. Another approach is to use a source term and an iterative procedure as proposed initially in [13]. In their work, they also derive a Chapman-Enskog analysis of the scheme. It was latter extended to phase change problems with convection, see [14, 15]. The same model is used in [16], with only one iteration for the source term at the expense of a lower accuracy. An implicit method is derived in [17], introducing a range of temperature for melting to model binary alloys which can be seen as a regularization of the enthalpy. Another approach consists in modifying the equilibrium function to recover enthalpy as proposed by [18]. Authors of [19] compare the hybrid LBM finite differences scheme introduced in [20] and the enthalpy approach in the context of phase change materials with convection. They found that the enthalpy method leads to more accurate results.

<footnote><sup>a</sup>Univ Rouen Normandie, CNRS, Normandie Univ, LMRS UMR 6085, F-76000 Rouen, France</footnote>

<footnote><sup>b</sup>Univ Angers, CNRS, LAREMA - UMR 6093, SFR MATHSTIC, F-49000 Angers, France</footnote>

<footnote>*Corresponding author.</footnote>

<footnote>Email addresses: francky.luddens@univ-rouen.fr (Francky Luddens), corentin.lothode@univ-angers.fr (Corentin Lothodé), ionut.danaila@univ-rouen.fr (Ionut Danaila)</footnote>
As noted in [21], the enthalpy method generates two non-physical effects when dealing with pure materials with a single fusion/solidification temperature: (i) the phase change interface increases step by step like a staircase; (ii) the temperature history at a given spatial point leads to a number of rogue plateaux in addition to the physical temperature plateau.

Following the recent work of [22], we propose another formulation of the source term, mixing a local Newton algorithm and a regularized enthalpy method. This approach has already shown promising results in the context of numerical methods based on finite elements. Adapting this method to LBM while keeping efficient parallelization is not a trivial task since it involves solving a non-linear problem during the stream-collide step.

The paper is organized as follows: in §2 we describe the Stefan problem. §3 provides a description of the discretization, and the link between a mesh (e.g. suitable for finite differences schemes) and the lattice is expressed. We also present previous methods to handle the enthalpy term. The novelty of this paper is described in §4, with a Newton method to compute the regularized enthalpy as a source term. §5 is devoted to the numerical analysis of the method, using the Chapman-Enskog expansion. Finally, §6 contains numerical illustrations in 1D and 2D to assess the stability and accuracy of the method.

## 2. Stefan problem and total enthalpy model

Consider the melting of a semi-infinite one-dimensional pure material, initially at uniform temperature $T_0 < T_f$, where $T_f$ is the fusion temperature. At the left border $x = 0$, the temperature is suddenly increased to $T_h > T_f$ and maintained constant for any time instant $t > 0$. The liquid phase starts to propagate from left to right. Solid and liquid domains are denoted by subscripts $s$ and $l$, respectively. They are separated by the interface located at $x_f(t)$. We use the following scaling:

* $x \to x/L$,

* $t \to t/(L^2/\alpha_l)$,

* $\theta = (T - T_f)/(T_h - T_f)$,

with $L$ a reference length scale and $\alpha_l = k_l/(\rho_l c_l)$ the thermal diffusivity. Density $\rho$, specific heat $c$ and thermal conductivity $k$ are constants, but different from one phase to another.

The Stefan problem is completely defined by writing the corresponding heat equation in each domain, with appropriate boundary conditions, and by imposing the balance of heat fluxes at the interface (Stefan condition). In the non-dimensional setting, the problem formulation is summarized in Eq. (1), where $h_{sl}$ is the latent heat and $Ste$ the Stefan number.

$$ \left. \begin{array}{l} \text{Liquid phase} \\ \frac{\partial \theta_l}{\partial t} = \frac{\partial^2 \theta_l}{\partial x^2} \\ \theta_l(0, t) = 1, \quad \theta_l(x_f, t) = 0 \end{array} \right| \quad \text{Solid phase} \quad \begin{array}{l} \frac{\partial \theta_s}{\partial t} = \alpha \frac{\partial^2 \theta_s}{\partial x^2}, \quad \alpha = \frac{\alpha_s}{\alpha_l}, \\ \theta_s(x_f, t) = 0, \quad \theta_s(+\infty, t) = \theta_0 < 0 \end{array} \qquad (1) $$

Stefan condition at the interface $x = x_f$

$$ -\frac{\partial \theta_l}{\partial x} + k \frac{\partial \theta_s}{\partial x} = \frac{\rho}{Ste} \frac{dx_f}{dt}, \quad \text{with} \quad \rho = \frac{\rho_s}{\rho_l}, \quad k = \frac{k_s}{k_l}, \quad Ste = \frac{c_l(T_h - T_f)}{h_{sl}}. $$

The following solution to Eq. (1) is obtained by assuming that the interface is moving following the law $x_f = 2\lambda\sqrt{t}$, with $\lambda$ a constant:

$$ \theta_l(x, t) = 1 - \frac{1}{\operatorname{erf}(\lambda)} \operatorname{erf}\left(\frac{x}{2\sqrt{t}}\right), \quad \theta_s(x, t) = \theta_0 - \frac{\theta_0}{\operatorname{erfc}\left(\frac{\lambda}{\sqrt{\alpha}}\right)} \operatorname{erfc}\left(\frac{x}{2\sqrt{\alpha t}}\right), \qquad (2) $$
where $erf(z) = (2/\sqrt{\pi}) \int_{0}^{z} \exp(-\tau^2)d\tau$ is the error function and $erfc(z) = 1 - erf(z)$. The constant $\lambda$ is determined from the Stefan condition, i. e. by solving the following equation:

$$ \frac{\exp(-\lambda^2)}{\operatorname{erf}(\lambda)} + \theta_0 \frac{k \exp(-\frac{\lambda^2}{\alpha})}{\sqrt{\alpha} \operatorname{erfc}(\frac{\lambda}{\sqrt{\alpha}})} = \frac{\rho \lambda \sqrt{\pi}}{\operatorname{Ste}}. \quad (3) $$

The total enthalpy model used to solve numerically the Stefan problem (1) is based on a temperature transformation of the energy conservation law. For the sake of clarity, we only present the case where the materials properties are similar in the liquid phase and the solid phase (hence $\rho = k = \alpha = 1$). The enthalpy $H$ is written as the sum of the sensible heat and the latent heat:

$$ H(\theta) = \theta + \frac{\varphi(\theta)}{\operatorname{Ste}} = \theta + \begin{cases} 0 & \text{if } \theta \le 0, \\ \frac{1}{\operatorname{Ste}} & \text{if } \theta > 0, \end{cases} $$

with $\varphi$ the liquid fraction (taking the value 1 if $\theta > 0$ and 0 if $\theta \le 0$). The following equation is obtained

$$ \frac{\partial H}{\partial t} - \Delta \theta = 0, \quad (4) $$

which is equivalent to

$$ \frac{\partial \theta}{\partial t} - \Delta \theta + \frac{1}{\operatorname{Ste}} \frac{\partial \varphi(\theta)}{\partial t} = 0. \quad (5) $$

Note that the formulations (4) or (5) of the Stefan problem are valid in 2D and 3D for phase change problems without convection.

## 3. LBM methods for the Stefan problem

### 3.1. Mesh and lattice

In this note, we will use the LBM to solve the PDE (5). When needed, we will link the lattice to an actual mesh of the domain $D$, resulting in a uniform Cartesian mesh. We denote by $N$ the size of the lattice, $\delta x = \frac{1}{N}$, and $\delta t$ the time step. We define the lattice $\mathcal{L} = [0, N]^d$, where $d$ is the space dimension. To properly perform the Chapman-Enskog analysis, we introduce the following rescaling:

$$ t_L := \frac{t}{\delta t}, \quad \mathbf{x}_L := \frac{\mathbf{x}}{\delta x}, \quad \theta_L(t_L, \mathbf{x}_L) := \theta(t, \mathbf{x}), \quad H_L(t_L, \mathbf{x}_L) := H(t, \mathbf{x}). $$

Note that if $\mathbf{x}$ is a node on the mesh, $\mathbf{x}_L$ is a node on the lattice. For derivatives, the subscript $L$ means spatial differentiation with respect to $\mathbf{x}_L$. In this setting, Eqs. (4) and (5) read

$$ \partial_{t_L} H_L - \frac{\delta t}{\delta x^2} \Delta_L \theta_L = 0, \quad (6) $$

$$ \partial_{t_L} \theta_L - \frac{\delta t}{\delta x^2} \Delta_L \theta_L = - \frac{1}{\operatorname{Ste}} \partial_{t_L} \varphi(\theta_L). \quad (7) $$

### 3.2. LBM setting

Let us denote by $e_i$, $i = 0, \dots, q-1$ the set of discrete velocities for the lattice we want to consider. We denote by $w_i$ the standard weights associated to the lattice. We introduce the distribution functions $f_i$, $f_i^{eq}$, $i = 0, \dots, q$; $f_i^{eq}$ is the equilibrium distribution and its value depends on the selected scheme. The lattices associated to a $d$-dimensional problem with $q$ discrete velocities is denoted by $DdQq$. In particular, we will use D1Q3 lattices in 1D and D2Q5 in 2D.
Throughout this note, if $S$ is a function of $t_L$ and $\mathbf{x}_L$, or if $S_i$ is a distribution function, we set:

$$
\begin{align*}
S &= S(t_L, \mathbf{x}_L), \\
\tilde{S} &= S(t_L + 1, \mathbf{x}_L), \\
S_i^* &= S(t_L + 1, \mathbf{x}_L + \mathbf{e}_i), \\
\hat{S}_i &= S(t_L, \mathbf{x}_L + \mathbf{e}_i), \\
\tilde{S}_i &= S(t_L, \mathbf{x}_L - \mathbf{e}_i).
\end{align*}
$$

We denote by $g_i$ a source term (to be specified later). Following the Bhatnagar-Gross-Krook (BGK) approximation, the LBM algorithm reads:

$$
f_i^* - f_i = -\frac{1}{\tau} (f_i - f_i^{eq}) + g_i, \quad (8)
$$

where $\tau$ is a relaxation parameter.

For a distribution function $s_i$, we denote by $M_k(s)$ the moment of order $k$, i.e.

$$
\mathcal{M}_k(s) := \sum_{i=0}^{q-1} \underbrace{\mathbf{e}_i \otimes \cdots \otimes \mathbf{e}_i}_{k \text{ times}} s_i.
$$

To alleviate the notations, we also define:

$$
\Pi_k := \mathcal{M}_k(f), \quad \Lambda_k := \mathcal{M}_k(f^{eq}), \quad \mathcal{G}_k := \mathcal{M}_k(g).
$$

3.3. Implicit liquid fraction based method (ILFBM)

One of the first specific LB models for the total enthalpy equation (5) was introduced in [13]. It relies on the use of a source/sink term. The distribution $f_i$ is defined such that $\theta_L = M_0(f)$ and $f_i^{eq} = w_i \theta_L$. The source term $g_i$ is intended to satisfy

$$
g_i = \frac{w_i}{Ste} (\ell - \tilde{\ell}),
$$

where the additional field $\ell$ (resp. $\tilde{\ell}$) is an approximation of the liquid fraction at time $t$ (resp. $t + \delta t$). The goal is to recover $\mathcal{G}_0 = -\frac{1}{Ste}\partial_t\varphi(\theta_L)$, at least, up to third order in $\varepsilon$ (which is the parameter used in the Chapman-Enskog expansion). $\ell$ is computed from the enthalpy $H$ using:

$$
\ell = l(H) = \begin{cases} 0 & \text{if } H < 0, \\ \text{Ste} H & \text{if } 0 \le H \le \frac{1}{\text{Ste}}, \\ 1 & \text{if } H > \frac{1}{\text{Ste}}. \end{cases}
$$

The time step procedure to get $\tilde{f}_i$ and $\tilde{\ell}$ from $f_i$ and $\ell$ is given in Algorithm 1. The target equation (7) can be recovered from Chapman-Enskog expansion, see [13]. We refer to this approach as the Implicit Liquid Fraction Based Method (ILFBM).

3.4. Explicit enthalpy based method (EEBM)

Another approach was introduced in [18] to take into account the enthalpy and avoid iterations. $f_i^{eq}$ is modified
such that its zeroth order moment is exactly the enthalpy $H$; the other moments are unchanged. $\theta$ can be recovered
from $H$ using the following expression:

$$
\theta = \Theta(H) = \begin{cases} H, & \text{if } H < 0, \\ 0, & \text{if } 0 \le H \le \frac{1}{\text{Ste}}, \\ H - \frac{1}{\text{Ste}}, & \text{if } H > \frac{1}{\text{Ste}}, \end{cases}
$$

hence no source term is required. The time step procedure to get $\tilde{f}_i$ and $\tilde{\ell}$ from $f_i$ and $\ell$ is given in Algorithm 2. No inner loop is required, since $\theta$ can be explicitly computed from $H$. The target equation (6) can be recovered from Chapman-Enskog expansion, see [18]. We refer to this approach as the Explicit Enthalpy Based Method (EEBM).
**Algorithm 1** Implicit Liquid Fraction Based Method step

1: moments: $\theta \leftarrow \mathcal{M}_0(f)$ and $H \leftarrow \varphi(\theta)$

2: equilibrium: $f_i^{eq} \leftarrow w_i \theta$

3: inner loop initialization: $\ell^{(0)} \leftarrow \ell$, $f_i^{(0)} \leftarrow f_i$

4: **while** ($\ell^{(k)}$ and $f_i^{(k)}$ not converged) **do**

5: collision: $\Omega_i \leftarrow f_i - \frac{1}{\tau}(f_i - f_i^{eq}) + \frac{w_i}{Ste}(\ell - \ell^{(k-1)})$

6: streaming: $f_i^{(k)} \leftarrow \tilde{\Omega}_i$

7: update enthalpy: $H^{(k)} \leftarrow \mathcal{M}_0(f^{(k)}) + \frac{1}{Ste} \ell^{(k-1)}$

8: update liquid fraction : $\ell^{(k)} \leftarrow l(H^{(k)})$

9: **end while**

10: output: $\tilde{f}_i \leftarrow f_i^{(k)}$ and $\tilde{\ell} = \ell^{(k)}$

**Algorithm 2** Explicit Enthalpy Based Method step

1: moments: $H \leftarrow \mathcal{M}_0(f)$, $\theta \leftarrow \Theta(H)$

2: **if** ($i \neq 0$) **then**

3: equilibrium: $f_i^{eq} = w_i \theta$

4: **else**

5: equilibrium: $f_0^{eq} = H - (1 - w_0)\theta$

6: **end if**

7: collision: $\Omega_i = f_i - \frac{1}{\tau}(f_i - f_i^{eq})$

8: streaming: $\tilde{f}_i \leftarrow \tilde{\Omega}_i$

## 4. Implicit regularized enthalpy based method (IREBM)

### 4.1. Main idea of the scheme

In this section, we present our approach to compute the enthalpy, which we will refer to as the Implicit Regularized Enthalpy Based Method (IREBM). We first replace the discontinuous liquid fraction $\varphi$ in (7) with the regularized function

$$ \varphi_{\delta}(\theta) = \frac{1}{2} \left( 1 + \tanh\left(\frac{\theta}{\delta}\right) \right), $$

where $\delta$ is a small positive smoothing parameter (see [22]). We derive a LB model to recover the regularized equation

$$ \partial_{t_L} \theta_L - \frac{\delta t}{\delta x^2} \Delta_L \theta_L = - \frac{1}{Ste} \partial_{t_L} \varphi_{\delta}(\theta_L). \quad (9) $$

As for the ILFBM, we will use a source term and we want that $\theta_L = \mathcal{M}_0(f)$ and $f_i^{eq} = w_i \theta_L$. The idea of the algorithm is to execute an inner loop in order to set

$$ g_i = \hat{s}_i - s_i^*, \quad \text{with} \quad s_i = \frac{w_i}{Ste} \varphi_{\delta}(\theta_L). \quad (10) $$

Note that both explicit ($\hat{s}_i$) and implicit ($s_i^*$) terms are evaluated at $x_L + e_i$. The reason is that the use of $s_i^*$ instead of $\tilde{s}_i$ allows us to recover a local (in space) Newton problem. Then the use of $\hat{s}_i$ is required instead of a more usual $s_i$ in order to avoid extra terms in the Chapman-Enskog expansion, see (21).

The LBM update scheme is

$$ f_i^* - f_i = -\frac{1}{\tau}(f_i - f_i^{eq}) + (\hat{s}_i - s_i^*). \quad (11) $$

The only unknown in the right hand side is $s_i^*$, for which we would need $\tilde{\theta}_L$. Let us introduce $\Omega_i$ and $q_i$ such that

$$ \hat{q}_i = \Omega_i := \left(1 - \frac{1}{\tau}\right) f_i + \frac{1}{\tau} f_i^{eq}. \quad (12) $$
Then the previous equation reads

$$f_i(t_L + 1, \mathbf{x}_L) = q_i(t_L, \mathbf{x}_L) + \frac{w_i}{Ste} (\varphi_\delta(\theta(t_L, \mathbf{x}_L)) - \varphi_\delta(\theta(t_L + 1, \mathbf{x}_L))). \quad (13)$$

Taking the zeroth order moment of (13) yields:

$$\tilde{\theta}_L = \mathcal{M}_0(q) + \frac{1}{Ste} (\varphi_\delta(\theta_L) - \varphi_\delta(\tilde{\theta}_L)). \quad (14)$$

We solve this (local in $\mathbf{x}_L$) non linear equation using a Newton algorithm in order to get $\tilde{\theta}_L$. The complete algorithm to get $\tilde{f}_i$ from $f_i$ is given in Algorithm 3.

**Algorithm 3 Implicit Regularized Enthalpy Based Method step**

1: moments: $\theta \leftarrow \mathcal{M}_0(f)$

2: equilibrium: $f_i^{eq} \leftarrow w_i \theta$

3: collision: $\Omega_i \leftarrow (1 - \frac{1}{\tau}) f_i + \frac{1}{\tau} f_i^{eq}$

4: streaming: $q_i \leftarrow \tilde{\Omega}_i$

5: Newton initialization: $m_0 \leftarrow \mathcal{M}_0(q)$, $\theta^{(0)} = \theta$

6: **while** ($\theta^{(k)}$ not converged) **do**

7: $\theta^{(k+1)} \leftarrow \frac{m_0 + \frac{1}{Ste}\varphi_\delta(\theta) + \frac{1}{Ste}\varphi'_\delta(\theta^{(k)})\theta^{(k)} - \frac{1}{Ste}\varphi_\delta(\theta^{(k)})}{1 + \frac{1}{Ste}\varphi'_\delta(\theta^{(k)})}$

8: **end while**

9: source term: $g_i \leftarrow \frac{w_i}{Ste}(\varphi_\delta(\theta) - \varphi_\delta(\theta^{(k)}))$

10: update: $\tilde{f}_i = q_i + g_i$

## 4.2. Treatment of boundary conditions

In order to get a complete LBM scheme in 1D or 2D, one has to impose boundary conditions on the distribution $f_i$. For EEBM and IFLBM, this can be done using standard techniques for Dirichlet or Neumann boundary conditions, cf. e. g. [23, 24]. For IREBM, owing to (12) and (13), we have to set boundary condition on the distribution $q_i$. One can notice that $q_i(t_L, \mathbf{x}_L)$ is not defined if $\mathbf{x}_L - \mathbf{e}_i$ is not a point of the lattice, but this value is still required to get the updated distribution function. Since the algorithm depends on the lattice used to solve the numerical problem, we present some examples for simple lattices in 1D and 2D. Denote by $\Gamma_D$ (resp. $\Gamma_N$) the part of the boundary on which we want to impose the Dirichlet condition $\theta = \theta_{Dir}$ (resp. the Neumann condition $\partial_n \theta = 0$).

### 4.2.1. Dirichlet boundary condition on a D1Q3 or D2Q5 lattice

Let us consider $\mathbf{x}_L \in \Gamma_D$ such that $\mathbf{x}_L$ is not on a corner (the idea can be adapted for the corners). From (14) and the definition of the moments, we obtain, for this point $\mathbf{x}_L$, that

$$\tilde{\theta} = \sum_{i=0}^{q-1} q_i + \frac{1}{Ste} (\varphi_\delta(\theta) - \varphi_\delta(\tilde{\theta})).$$

Using the Dirichlet boundary condition, this can be rewritten as:

$$\mathcal{M}_0(q) = \tilde{\theta}_{Dir} + \frac{1}{Ste} (\varphi_{\delta}(\theta_{Dir}) - \varphi_{\delta}(\tilde{\theta}_{Dir})). \quad (15)$$

For a general lattice, one might use any technique to impose a Dirichlet condition on $q_i$, according to (15). In the case of a D1Q3 or D2Q5 lattice, only one of the $q_i(t, \mathbf{x}_L)$ is not well-defined in (12). Let us denote by $i_x$ its index. Then (15) reads

$$q_{i_x} = \tilde{\theta}_{Dir} + \frac{1}{Ste} (\varphi_{\delta}(\theta_{Dir}) - \varphi_{\delta}(\tilde{\theta}_{Dir})) - \sum_{i \neq i_x} q_i. \quad (16)$$

Since all the terms in the right hand side are well defined and known, this is a consistent way of defining $q_{i_x}$.
### 4.2.2. Bounce-back condition

A standard way to impose the homogeneous Neumann boundary condition is to use the so-called "bounce-back" condition [25]. For any $\mathbf{x}_L \in \Gamma_N$, we denote by $i_x$ any index such that $q_{i_x}(t, \mathbf{x}_L)$ is not well-defined in (12). The bounce-back condition states that

$$f_{i_x} = f_{\bar{i}_x},$$

at any time, where $\bar{i}_x$ is the index of the opposite direction to the discrete velocity $e_{i_x}$. Using (13), this leads to

$$q_{i_x} + \frac{w_{i_x}}{Ste} (\varphi_{\delta}(\theta) - \varphi_{\delta}(\tilde{\theta})) = q_{\bar{i}_x} + \frac{w_{\bar{i}_x}}{Ste} (\varphi_{\delta}(\theta) - \varphi_{\delta}(\tilde{\theta})).$$

Since $w_{i_x} = w_{\bar{i}_x}$, this reduces to a bounce-back condition on $q$, i.e. $q_{i_x} = q_{\bar{i}_x}$. Since the term in the right-hand side is well defined and known, this is a consistent way of defining $q_{i_x}$. Note also that this bounce-back condition can be used on any lattice.

### 4.2.3. Neumann boundary condition on a D2Q5 lattice

Let us consider $\mathbf{x}_L \in \Gamma_N$ such that $\mathbf{x}_L$ is not a corner (the idea can be adapted for the corners). In the particular case of a D2Q5 lattice, only one of the $q_i(t, \mathbf{x}_L)$ is not well-defined in (12). Let us denote by $i_x$ its index. An alternative solution to the bounce back condition (with order 1 accuracy in space) is to impose (with our notations) that $\theta(t_L, \mathbf{x}_L) = \theta(t_L, \mathbf{x}_L + \mathbf{e}_{i_x})$ for any time. From (14), we have

$$\tilde{\theta}(\mathbf{x}_L) = \sum_{i=0}^{q-1} q_i(\mathbf{x}_L) + \frac{1}{Ste} (\varphi_{\delta}(\theta(\mathbf{x}_L)) - \varphi_{\delta}(\tilde{\theta}(\mathbf{x}_L))),$$

$$\tilde{\theta}(\mathbf{x}_L + \mathbf{e}_{i_x}) = \sum_{i=0}^{q-1} q_i(\mathbf{x}_L + \mathbf{e}_{i_x}) + \frac{1}{Ste} (\varphi_{\delta}(\theta(\mathbf{x}_L + \mathbf{e}_{i_x})) - \varphi_{\delta}(\tilde{\theta}(\mathbf{x}_L + \mathbf{e}_{i_x}))).$$

The applied Neumann condition at time $t$ ensures that $\theta(\mathbf{x}_L) = \theta(\mathbf{x}_L + \mathbf{e}_{i_x})$. Thus, imposing $\tilde{\theta}(\mathbf{x}_L) = \tilde{\theta}(\mathbf{x}_L + \mathbf{e}_{i_x})$ is equivalent to satisfy

$$\sum_{i=0}^{q-1} q_i(\mathbf{x}_L) = \sum_{i=0}^{q-1} q_i(\mathbf{x}_L + \mathbf{e}_{i_x}).$$

As a result, we have

$$q_{i_x}(\mathbf{x}_L) = \sum_{i=0}^{q-1} q_i(\mathbf{x}_L + \mathbf{e}_{i_x}) - \sum_{i \neq i_x} q_i(\mathbf{x}_L). \qquad (17)$$

Since all the terms in the right-hand side are well defined and known, this is a consistent way of defining $q_{i_x}(\mathbf{x}_L)$.

## 5. Numerical analysis of the scheme

The goal of this section is to recover the target equation (7) from the LBM setting, using a Chapman-Enskog expansion. We start from the LBM update scheme (11) and assume that the following holds (standard assumptions for the Chapman-Enskog expansion), for a small parameter $\epsilon$:

$$f_i = f_i^{(0)} + \epsilon f_i^{(1)} + \epsilon^2 f_i^{(2)},$$

$$f_i^{eq} = f_i^{eq,0},$$

$$\partial_{t_L} = \epsilon \partial_{t_1} + \epsilon^2 \partial_{t_2},$$

$$\nabla_L = \epsilon \nabla_{\bar{x}}.$$

Since we use that $f_i^{eq} = w_i \theta_L$, we obtain

$$\Lambda_0 = \Pi_0 = \theta_L, \qquad (18)$$

$$\Lambda_1 = 0, \qquad (19)$$

$$\Lambda_2 = \frac{1}{3} \theta_L \mathbb{I}. \qquad (20)$$
We also denote by $\Pi_k^{(m)}$ the k-th order moment of the distribution $f_i^{(m)}$. Expanding (11) up to third order in $\epsilon$, we get:

$$ \epsilon D_{1,i} f_i + \frac{\epsilon^2}{2} D_{1,i}^2 f_i + \partial_{t_2} f_i = -\frac{1}{\tau} (f_i - f_i^{eq}) - \left( \epsilon \partial_{t_1} s_i + \epsilon^2 \partial_{t_2} s_i + \frac{\epsilon^2}{2} \partial_{t_1}^2 s_i + \epsilon^2 \partial_{t_1} (e_i \cdot \nabla_{\bar{x}} s_i) \right), \quad (21) $$

where $D_{1,i}f = \partial_{t_1}f + e_i \cdot \nabla_{\bar{x}}f$. By identification with respect to consecutive orders of $\epsilon$, we obtain:

$$ \epsilon^0 : f_i^{(0)} = f_i^{eq}, \quad (22) $$

$$ \epsilon^1 : D_{1,i} f_i^{(0)} = -\frac{1}{\tau} f_i^{(1)} - \partial_{t_1} s_i, \quad (23) $$

$$ \epsilon^2 : \partial_{t_2} f_i^{(0)} + D_{1,i} f_i^{(1)} + \frac{1}{2} D_{1,i}^2 f_i^{(0)} = -\frac{1}{\tau} f_i^{(1)} - \left( \partial_{t_2} s_i + \frac{1}{2} \partial_{t_1}^2 s_i + \partial_{t_1} (e_i \cdot \nabla_{\bar{x}} s_i) \right). \quad (24) $$

Taking the moments of (22) yields $\Pi_k^{(0)} = \Lambda_k$, for any $k$. Since we have $\Pi_0 = \Lambda_0 = \theta_L$, we infer that, for any $l \ge 1$:

$$ \Pi_0^{(l)} = 0. \quad (25) $$

Injecting (23) in (24), we obtain:

$$ \partial_{t_2} f_i^{(0)} + \left(1 - \frac{1}{2\tau}\right) D_{1,i} f_i^{(1)} = -\frac{1}{\tau} f_i^{(1)} - \left(\partial_{t_2} s_i + \frac{1}{2} \partial_{t_1}^2 s_i + \partial_{t_1} (e_i \cdot \nabla_{\bar{x}} s_i)\right) + \frac{1}{2} D_{1,i} \partial_{t_1} s_i. \quad (26) $$

Taking the zeroth and the first order moment of (23) and the zeroth order moment of (26) leads to:

$$ \partial_{t_1} \theta_L = - \frac{1}{Ste} \partial_{t_1} \varphi_\delta(\theta_L), \quad (27) $$

$$ \frac{1}{3} \nabla_{\bar{x}} \cdot (\theta_L \mathbb{I}) = - \frac{1}{\tau} \Pi_1^{(1)}, \quad (28) $$

$$ \partial_{t_2} \theta_L + \left(1 - \frac{1}{2\tau}\right) \nabla_{\bar{x}} \cdot \Pi_1^{(1)} = - \frac{1}{Ste} \partial_{t_2} \varphi_\delta(\theta_L). \quad (29) $$

Inserting (28) in (29) gives:

$$ \partial_{t_2} \theta_L - \frac{1}{3} \left( \tau - \frac{1}{2} \right) \nabla_{\bar{x}}^2 \theta_L = - \frac{1}{Ste} \partial_{t_2} \varphi_\delta(\theta_L). \quad (30) $$

We multiply (27) and (30) by $\epsilon$ and $\epsilon^2$ respectively and sum these equations to get:

$$ \partial_{t_L} \theta_L - \frac{1}{3} \left( \tau - \frac{1}{2} \right) \Delta_L \theta_L = - \frac{1}{Ste} \partial_{t_L} \varphi_\delta(\theta_L), \quad (31) $$

up to the third order in $\epsilon$. We can see that (31) corresponds to (7), provided that:

$$ \tau - \frac{1}{2} = \frac{3\delta t}{\delta x^2}. $$

# 6. Numerical illustrations

## 6.1. 1D problem

To validate the algorithm developed in §4, namely the implicit regularized enthalpy based method (IREBM), we consider a case from [18]. We compare results obtained by two other algorithms, the explicit enthalpy based method (EEBM) from [18] and the Implicit Liquid Fraction Based Method (ILFBM) from [13]. We use the numerical methods described in §2 and the exact solution as reference. We consider a one dimensional bar of length L made of aluminum. The temperature is fixed at both ends, with $T_h$ on the left and $T_c$ on the right. During the simulation, physical quantities are rendered dimensionless, and dimensionless temperature is denoted by $\theta$. The case is initialized with a temperature
distribution satisfying the Stefan condition for the case where most of the bar is solid. More precisely, the position of the interface is given by $x_f(t) = 2\lambda\sqrt{t}$, where $\lambda$ is computed from (3). We choose $t_0$ such that $x_f(t = t_0) = 0.01$. Starting from the initial temperature distribution given by (2), we track the displacement of the solid-liquid interface. One can refer to Figure 1 for a sketch of the simulation; physical properties are given in Table 1.

<figure><img src="image_3.png" /><figcaption>Figure 1: Sketch of the simulation configuration.</figcaption></figure>

<table><caption>Table 1: Physical quantities and values used for the 1D case.</caption><thead><tr><th>physical quantity</th><th>symbol</th><th>value</th><th>unit</th></tr></thead><tbody><tr><td>density</td><td>ρ<sub>l</sub> = ρ<sub>s</sub></td><td>2500</td><td>kg·m<sup>-3</sup></td></tr><tr><td>specific heat capacity</td><td>c<sub>l</sub> = c<sub>s</sub></td><td>1</td><td>kJ·kg<sup>-1</sup>·K<sup>-1</sup></td></tr><tr><td>thermal conductivity</td><td>k<sub>l</sub> = k<sub>s</sub></td><td>0.2</td><td>kW·m<sup>-1</sup>·K<sup>-1</sup></td></tr><tr><td>thermal diffusivity</td><td>α<sub>l</sub> = α<sub>s</sub></td><td>8 × 10<sup>-5</sup></td><td>m<sup>2</sup>·s<sup>-1</sup></td></tr><tr><td>latent heat</td><td>h<sub>sl</sub></td><td>350</td><td>kJ·kg<sup>-1</sup>·K<sup>-1</sup></td></tr><tr><td>fusion temperature</td><td>T<sub>f</sub></td><td>500</td><td>°C</td></tr><tr><td>temperature left side</td><td>T<sub>h</sub></td><td>600</td><td>°C</td></tr><tr><td>temperature right side</td><td>T<sub>c</sub></td><td>450</td><td>°C</td></tr><tr><td>reference length</td><td>L</td><td>1</td><td>m</td></tr></tbody></table>

The non-dimensional formulation of this problem is obtained by using the as reference temperature the difference $\delta T = T_h - T_f$, resulting in the Stefan number $Ste = 0.2857$. The initial temperature profile is shown in Fig. 2 and exhibits a discontinuity of derivatives at $x_f$, as expected since the phase change occurs at this location.

<figure><img src="image_4.png" /><figcaption>Figure 2: Initial solution for the 1D case.</figcaption></figure>

In Fig. 3, we compare the three methods using the same spatial and time resolutions. EEBM and ILFBM behave almost identically, resulting in an oscillation of the error on the position of the interface between $2 \times 10^{-3}$ and $1 \times 10^{-2}$. For the implicit regularized enthalpy based method (IREBM), the error on the position of the interface oscillates
between $3 \times 10^{-3}$ and $4 \times 10^{-3}$. Note that the error fluctuates a lot less with IREBM, which indeed is the result of the added regularization. One needs to keep in mind that the problem is modified when using the regularization, and methods EEBM, ILFBM will not converge to the solution of IREBM if $\delta$ is kept constant.

<figure><img src="image_5.png" /><figcaption>Figure 3: Position of the interface (left) and error with respect to the exact solution (right). Comparison between EEBM, ILFBM, IREBM (with $\delta = 0.005$) $N = 801$.</figcaption></figure>

In Fig. 4, we can see the effect of $\delta$ on the evolution of the position and the oscillations of error of the interface. Larger is $\delta$, smaller the oscillations of the error are. If $\delta$ is too large, the regularization will offset the solution to another problem, so $\delta$ needs to be tuned in order to keep the solution close to the original problem.

<figure><img src="image_6.png" /><figcaption>Figure 4: Position of the interface (left) and error with respect to the exact solution (right). Comparison of IREBM for three values of $\delta$ ($\delta_{min} = 0.005$), for $N = 801$.</figcaption></figure>
In Fig. 5, we can see the effect of the resolution on the solution, when keeping $\delta$ constant. We note that the error on the position of the interface diminishes as the resolution $N$ increases, and the oscillations are reduced.

<figure><img src="image_7.png" /><figcaption>Figure 5: Position of the interface (left) and error with respect to the exact solution (right). Comparison of different mesh resolutions for IREBM, with $\delta = 0.005$.</figcaption></figure>

In Fig. 6, we compare the three methods with the same spatial and time resolution. All methods exhibit very similar behaviors. As expected, the maximum of the error is near the phase change interface. The IREBM is having a slightly lower error than the two other methods. The maximum of the error is $4.06 \times 10^{-3}$ for the EEBM, $4.50 \times 10^{-3}$ for the ILFBM and $3.41 \times 10^{-3}$ for the IREBM.

<figure><img src="image_8.png" /><figcaption>Figure 6: Solution (left) and error (right) at final time ($t = 160s$) for all the methods (EEBM, ILFBM, IREBM with $\delta = 0.005$) and $N = 801$.</figcaption></figure>

Finally, we report the impact of Newton iterations on the computational time in Tab. 2. As expected, IREBM is
slower than EEBM because of Newton iterations. However, IREBM is significantly faster than ILBM for two reasons: (i) the Newton algorithm usually requires a low number of iterations to converge; (ii) inner iterations are performed on only one field, while the ILFBM inner loop involves all *q* fields of the particle distribution function and a streaming process. For these reasons, the extra cost of the Newton loop is expected to be even lower in 2D or 3D, when larger values of *q* in the lattice are necessary.

<table><caption>Table 2: 1D case, lattice D1Q3: computational time (in seconds, sequential code on Mac M1), with $\tau = 0.62$ and $\delta = 0.01$ for IREBM.</caption><thead><tr><th>N</th><th>EEBM</th><th>IREBM</th><th>ILFBM</th></tr></thead><tbody><tr><td>201</td><td>0.026</td><td>0.106</td><td>0.422</td></tr><tr><td>401</td><td>0.200</td><td>0.739</td><td>2.997</td></tr><tr><td>801</td><td>1.538</td><td>5.586</td><td>22.01</td></tr><tr><td>1601</td><td>12.43</td><td>43.79</td><td>164.06</td></tr></tbody></table>

## 6.2. 2D problem

In order to show the 2D capabilities of the developed algorithm, we compare our results to [26]. They based their comparison on an approximate *analytical* solution from [27]. Other analysis of solutions at corners is given in [28]. The setup consists of a squared cavity of $8m$ wide, initially at $T = T_i = 10°C$ and cooled at all sides to $T = T_Γ = -20°C$. Using the symmetries of the problem, we solve Eq. 5 in a square representing the lower-left quarter of the original domain, see Fig. 8. One can refer to Fig. 7 for a sketch of the simulation; physical properties are given in Tab. 3.

The computational domain (see e.g. Figure 8) is then a square, in which Eq. (5) is solved. The bottom and left walls are kept at the cold temperature, resulting in a Dirichlet boundary condition. To take into account symmetry properties, the system is supplemented with the homogeneous Neumann condition on the top and right walls.

<figure><img src="image_9.png" /><figcaption>Figure 7: Sketch of the 2D simulation configuration.</figcaption></figure>
<table><caption>Table 3: Physical quantities and their values used for the 2D case.</caption><thead><tr><th>physical quantity</th><th>symbol</th><th>value</th><th>unit</th></tr></thead><tbody><tr><td>density</td><td>ρ<sub>l</sub> = ρ<sub>s</sub></td><td>1000</td><td>kg·m<sup>-3</sup></td></tr><tr><td>specific heat capacity</td><td>c<sub>l</sub> = c<sub>s</sub></td><td>1.762</td><td>kJ·kg<sup>-1</sup>·K<sup>-1</sup></td></tr><tr><td>thermal conductivity</td><td>k<sub>l</sub> = k<sub>s</sub></td><td>2.22×10<sup>-3</sup></td><td>kW·m<sup>-1</sup>·K<sup>-1</sup></td></tr><tr><td>thermal diffusivity</td><td>α<sub>l</sub> = α<sub>s</sub></td><td>1.26 × 10<sup>-6</sup></td><td>m<sup>2</sup>·s<sup>-1</sup></td></tr><tr><td>latent heat</td><td>h<sub>sl</sub></td><td>338</td><td>kJ·kg<sup>-1</sup>·K<sup>-1</sup></td></tr><tr><td>fusion temperature</td><td>T<sub>f</sub></td><td>0</td><td>°C</td></tr><tr><td>initial temperature</td><td>T<sub>i</sub></td><td>10</td><td>°C</td></tr><tr><td>temperature boundary</td><td>T<sub>Γ</sub></td><td>-20</td><td>°C</td></tr><tr><td>reference length</td><td>L</td><td>1</td><td>m</td></tr></tbody></table>

The non-dimensional formulation of this problem is obtained by using as reference temperature the difference $δT = T_i - T_f$, resulting in the Stefan number $Ste = 0.0521$. Figure 8 compares both EEBM and IREBM (ILFBM did not run successfully) at a final time $t = 20h$. Both methods behave similarly, excepting at the interface where EEBM, as expected, has an oscillatory behavior. The level set $T = -5°C$ is identical for both methods, but positive isovalues are more distinct.

<figure><img src="image_10.png" /><figcaption>Figure 8: 2D case: distribution of temperature for EEBM (left) and IREBM (center) and position of three isovalues (right). Simulations for $N = 201$ and $δ = 2 × 10^{-2} ≈ Δx$.</figcaption></figure>

To assess that the displacement of the isovalues is not due to the spatial resolution, Fig. 9 compares both EEBM and IREBM at a final time $t = 20h$ but for a higher resolution than previously ($n_{x,y} = 401$ compared to $n_{x,y} = 201$). The regularization value $δ$ is kept the same ($δ = 2 × 10^{-2}$). The same overall behavior is observed for both methods, with the same oscillating behavior for EEBM and $T = 0°C$.

<figure><img src="image_11.png" /><figcaption>Figure 9: 2D case: distribution of temperature for EEBM (left) and IREBM (center) and position of some isovalues (right). Simulations for $N = 401$ and $δ = 2 × 10^{-2} ≈ 2Δx$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10: 2D case: temperature along the *x* = *y* axis at *t* = 20h. Comparison of EEBM and IREBM for different mesh resolutions.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: 2D case: temperature along the *x* = *y* axis for *t* = 5, 20, 40h (IREBM simulation, *N* = 401).</figcaption></figure>

Figures 10 and 11 show the temperature along the *x* = *y* axis. In Fig. 10, EEBM and IREBM are compared for two discretizations (*N* = 201 and *N* = 401), and behave similarly. In Fig. 11, several time steps are shown. Our results match those of [26].

Finally, we report the impact of the Newton iterations on the computational time in Tab. 4. The IREBM is slower than the EEBM with a close to 3 for the D2Q5 lattice and 2.5 for the D2Q9 lattice. This is indeed better than the 1D case (see Tab. 2) where the ratio was close to 3.5. As expected, for large values of *q*, the extra cost for the Newton iterations is reduced. This is explained by the fact that the Newton iterations are performed on only one field, independently of the value of *q*. Note that ILFBM results are not shown here because the computational time is significantly larger. The reason is that, for all spatial resolutions and time steps considered here, the inner iterations failed to converge, and the scheme does not yield a proper solution.

Note also that, from a strict computational point of view, when *N* is multiplied by 2, the computational cost of a LBM scheme for a 2D case should be multiplied by 16. Here we notice that going from *N* = 201 to *N* = 401, the ratio is greater than 25. This is due to the fact that for *N* = 101 and *N* = 201, all the lattice stands in the L2 memory,
but for $N = 401$ this is no longer the case (12MB on Apple M1).

<table><caption>Table 4: 2D case: computational time (in seconds, sequential code on Mac M1), with $\tau = 0.84$ and $\delta = 0.01$ for IREBM.</caption><thead><tr><th rowspan="2">N</th><th colspan="3">D2Q5</th><th colspan="3">D2Q9</th></tr><tr><th>EEBM</th><th>IREBM</th><th>ratio</th><th>EEBM</th><th>IREBM</th><th>ratio</th></tr></thead><tbody><tr><td>101</td><td>0.064</td><td>0.250</td><td>3.90</td><td>0.104</td><td>0.303</td><td>2.91</td></tr><tr><td>201</td><td>0.875</td><td>3.794</td><td>4.34</td><td>1.697</td><td>5.133</td><td>3.02</td></tr><tr><td>401</td><td>25.194</td><td>74.245</td><td>2.95</td><td>43.452</td><td>106.58</td><td>2.45</td></tr><tr><td>801</td><td>656.58</td><td>1675.79</td><td>2.55</td><td>1028.56</td><td>2188.74</td><td>2.13</td></tr></tbody></table>

## 7. Conclusion

A new implicit regularized enthalpy based scheme was developed to solve the Stefan problem using the LBM formalism. The new method was analysed using both Taylor expansion and Chapman-Enskog analysis. This new approach (IREBM) incorporates ideas from both the implicit liquid fraction based method and also from the explicit enthalpy based method. The regularization introduced in the method reduces the oscillations of the position of the interface observed in the results obtained with previously published methods. Considering the position of the interface, the results obtained with IREBM show a better stability and lower error when compared to the two other methods. At a given time, the overall solution is similar to that obtained with EEBM or ILFBM. Parallelization and vectorization is kept identical and the arithmetic intensity is augmented due to the Newton loop. Future work will include new capabilities, new test cases on 3D configurations, and GPU open-source implementation.

## Acknowledgements

This project was co-financed by the European Union with the European regional development funds and by the Normandy Regional Council via the ELBA project. Part of this work was performed using computing resources of CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie, France). The authors also thank Pr. U. Rüde and his group at FAU (Friedrich-Alexander-Universität) for fruitful discussions and insights regarding the detailed implementation of LBM.

## References

[1] G. Lamé, B. P. Clapeyron, Memoire sur la solidification par refroidissement d'un globe liquide, Annales Chimie Physique 47 (1831) 250-256.

[2] D. A. Tarzia, A bibliography on moving-free boundary problems for the heat-diffusion equation. The Stefan and related problems., Materials 2 (2000) 1-297.

[3] A. Kowalewski, D. Gobin, Phase change with convection: modelling and validation, Springer, 2004.

[4] L. Rubinstein, The Stefan problem: Comments on its present state, IMA Journal of Applied Mathematics 24 (3) (1979) 259-277.

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

[6] C. Elliott, Error analysis of the enthalpy method for the Stefan problem, IMA journal of numerical analysis 7 (1) (1987) 61-71.

[7] C. H. Li, A finite-element front-tracking enthalpy method for Stefan problems, IMA Journal of Numerical Analysis 3 (1) (1983) 87-107.

[8] M. Muhieddine, E. Canot, R. March, Various approaches for solving problems in heat conduction with phase change, International Journal on Finite Volumes (2009) 19.

[9] S. Savović, J. Caldwell, Finite difference solution of one-dimensional Stefan problem with periodic boundary conditions, International journal of heat and mass transfer 46 (15) (2003) 2911-2916.

[10] D. Wolf-Gladrow, A lattice Boltzmann equation for diffusion, Journal of statistical physics 79 (5) (1995) 1023-1032.

[11] G. De Fabritiis, A. Mancini, D. Mansutti, S. Succi, Mesoscopic models of liquid/solid phase transitions, International Journal of Modern Physics C 9 (08) (1998) 1405-1415.

[12] W. Miller, S. Succi, D. Mansutti, Lattice Boltzmann model for anisotropic liquid-solid phase transition, Physical review letters 86 (16) (2001) 3578.

[13] W.-S. Jiaung, J.-R. Ho, C.-P. Kuo, Lattice Boltzmann method for the heat conduction problem with phase change, Numerical Heat Transfer: Part B: Fundamentals 39 (2) (2001) 167-187.

[14] D. Chatterjee, S. Chakraborty, An enthalpy-based lattice Boltzmann model for diffusion dominated solid-liquid phase transformation, Physics Letters A 341 (1-4) (2005) 320-330.
[15] D. Chatterjee, S. Chakraborty, An enthalpy-source based lattice Boltzmann model for conduction dominated phase change of pure substances, International Journal of Thermal Sciences 47 (5) (2008) 552–559.

[16] C. Huber, A. Parmigiani, B. Chopard, M. Manga, O. Bachmann, Lattice Boltzmann model for melting with natural convection, International Journal of Heat and Fluid Flow 29 (5) (2008) 1469–1480.

[17] M. Eshraghi, S. D. Felicelli, An implicit lattice Boltzmann model for heat conduction with phase change, International Journal of Heat and Mass Transfer 55 (9-10) (2012) 2420–2428.

[18] R. Huang, H. Wu, P. Cheng, A new lattice Boltzmann model for solid-liquid phase change, International Journal of Heat and Mass Transfer 59 (2013) 295–301.

[19] K. Luo, F.-J. Yao, H.-L. Yi, H.-P. Tan, Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials, Applied Thermal Engineering 86 (2015) 238–250.

[20] P. Lallemand, L.-S. Luo, Hybrid finite-difference thermal lattice Boltzmann equation, International Journal of Modern Physics B 17 (01n02) (2003) 41–47.

[21] V. Voller, M. Cross, Accurate solutions of moving boundary problems using the enthalpy method, International journal of heat and mass transfer 24 (3) (1981) 545–556.

[22] G. Sadaka, A. Rakotondrandisa, P.-H. Tournier, F. Luddens, C. Lothodé, I. Danaila, Parallel finite-element codes for the simulation of two-dimensional and three-dimensional solid–liquid phase-change systems with natural convection, Computer Physics Communications 257 (2020) 107492.

[23] G. Kaluza, The numerical solution of transient heat conduction problem using the lattice Boltzmann method, Sci. Res. Inst. Math. Comput. Sci. 1 (13) (2012) 23–30.

[24] Y. Chen, C. Müller, A Dirichlet boundary condition for the thermal lattice Boltzmann method, International Journal of Multiphase Flow 123 (2020) 103184.

[25] D. D'humières, P. Lallemand, Lattice gas automata for fluid mechanics, Physica A: Statistical Mechanics and its Applications 140 (1-2) (1986) 326–335.

[26] R. Prapainop, K. Maneeratana, Simulation of ice formation by the finite volume method, Simulation 26 (1) (2004) 56.

[27] L. Jiji, K. Rathjen, T. Drzewiecki, Two-dimensional solidification in a corner, International Journal of Heat and Mass Transfer 13 (1) (1970) 215–218.

[28] J. King, D. Riley, A. Wallman, Two–dimensional solidification in a corner, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 455 (1989) (1999) 3449–3470.