1. Introduction
In the past decades, the lattice Boltzmann (LB) method has been established as an effective alternative flow solver over the commonly used computational fluid dynamics (CFD) methods. The LB method translates macroscopic physics to mesoscopic interactions between hypothetical particles, where distribution functions present the probability of mass density in a mean field (Higuera, Succi & Benzi Reference Higuera, Succi and Benzi1989). Since it was developed as an improved version of the lattice gas automata method (McNamara & Zanetti Reference McNamara and Zanetti1988; Higuera & Jiménez Reference Higuera and Jiménez1989; Higuera & Succi Reference Higuera and Succi1989), the LB method has been rapidly proven to be a reliable computational tool for a wide range of physics (Benzi, Succi & Vergassola Reference Benzi, Succi and Vergassola1992; He, Chen & Doolen Reference He, Chen and Doolen1998a; Yu et al. Reference Yu, Mei, Luo and Shyy2003; Tang et al. Reference Tang, Li, Wang, He and Tao2006; Wang, Wang & Li Reference Wang, Wang and Li2006; Basati, Mohammadipour & Niazmand Reference Basati, Mohammadipour and Niazmand2021), with a bright future ahead (Succi Reference Succi2015).
The LB method is easy to code, can handle complex geometries and is very efficient for parallel computing. These features make the LB method an attractive choice for the simulation of complex flows. One of the problems that has gained the attention of numerous researchers is multiphase flow simulation. The colour method that Rothman & Keller (Reference Rothman and Keller1988) introduced was among the first LB models to simulate multiphase flows. The idea of particle interactions introduced in this model inspired much later research in this field. Shan & Chen (Reference Shan and Chen1993) developed an LB scheme using interparticle potentials to model multiphase and multicomponent fluid flows. Although this pseudo-potential model can simulate many phase separation phenomena (Li et al. Reference Li, Luo, Kang, He, Chen and Liu2016; Sayyari, Kharmiani & Esfahani Reference Sayyari, Kharmiani and Esfahani2019), it suffers from some drawbacks. As He & Doolen (Reference He and Doolen2002) showed in their work, the definition of interparticle potential in this model could lead to thermodynamic inconsistencies. Furthermore, the interfacial tension in the original model is fixed by choosing the strength of the interparticle potential and cannot be modified independently (Sbragaglia et al. Reference Sbragaglia, Benzi, Biferale, Succi, Sugiyama and Toschi2007). Swift, Osborn & Yeomans (Reference Swift, Osborn and Yeomans1995) pointed out the importance of thermodynamic consistency in multiphase LB simulations. They introduced a new LB model named the free energy method. In contrast to interparticle potential models, which do not have any mechanism to track the interface explicitly, the free energy method utilized a modified Cahn–Hilliard equation (Cahn & Hilliard Reference Cahn and Hilliard1958) to capture the interface in an explicit way. Using a new equation provides the freedom to control the surface tension. Continuous improvements have been made to the free energy LB model since its introduction in 1995. The original free energy model does not satisfy the Galilean invariance. Holdych et al. (Reference Holdych, Rovas, Georgiadis and Buckius1998) solved this problem by developing an improved free energy model. He, Chen & Zhang (Reference He, Chen and Zhang1999) proposed a kinetic-theory-based multiphase LB model to simulate multiphase flows at moderate density ratios. Inamuro et al. (Reference Inamuro, Ogata, Tajima and Konishi2004) and Lee & Lin (Reference Lee and Lin2005) extended the density ratio limitation. Zheng, Shu & Chew (Reference Zheng, Shu and Chew2006) improved the accuracy by exactly recovering the Cahn–Hilliard equation in LB formulation. An additional force was introduced by Li et al. (Reference Li, Luo, Gao and He2012) to prevent possible inconsistencies in the momentum equation at high velocities.
Modelling thermal flows with the LB method began with Alexander, Chen & Sterling (Reference Alexander, Chen and Sterling1993) and Qian (Reference Qian1993) in 1993. In the subsequent years, many thermal LB models have been developed and adapted to different thermal problems including one-component liquid–solid phase transitions. De Fabritiis et al. (Reference de Fabritiis, Mancini, Mansutti and Succi1998) introduced an LB model to describe flows with liquid–solid phase transition using a fictitious chemical reaction term. Jiaung, Ho & Kuo (Reference Jiaung, Ho and Kuo2001) described the liquid–solid phase transition by enthalpy formulation and used liquid fraction to determine the phase change interfacial position. A slightly modified version of this model was later used by Huber et al. (Reference Huber, Parmigiani, Chopard, Manga and Bachmann2008) to simulate pure-substance melting coupled with natural convection. Chakraborty & Chatterjee (Reference Chakraborty and Chatterjee2007) proposed a hybrid thermal LB model for studying convection–diffusion phenomena with liquid–solid phase change. Eshraghi & Felicelli (Reference Eshraghi and Felicelli2012) suggested treating the latent heat source implicitly by solving a group of linear equations in the collision step. By modifying the equilibrium distribution function for the temperature, Huang, Wu & Cheng (Reference Huang, Wu and Cheng2013) developed a new approach to treat the latent heat source term.
The liquid–solid transformation is also adapted in multicomponent simulations. The LB developed by Parmigiani et al. (Reference Parmigiani, Huber, Bachmann and Chopard2011) combines the multicomponent Shan–Chen method with the thermal LB model (Huber et al. Reference Huber, Parmigiani, Chopard, Manga and Bachmann2008) to simulate multiphase flows in porous media with melting and solidification. An enthalpy model based on percolation theory in porous media (Semma et al. Reference Semma, El Ganaoui, Bennacer and Mohamad2008) was used by Sun, Gong & Li (Reference Sun, Gong and Li2015) to simulate the solidification of a water droplet on a cold flat plate using the Shan–Chen method. In a cold environment with forced convection, Zhao et al. (Reference Zhao, Dong, Li and Dou2017) investigated the freezing of a free-falling droplet using an enthalpy-based LB model. Xiong & Cheng studied the impact and solidification of a droplet on a cold smooth (Reference Xiong and Cheng2018b) and rough (Reference Xiong and Cheng2018a) substrate. The frost crystal nucleation process involved in liquid–solid phase change on a cold surface is studied by Gong et al. (Reference Gong, Hou, Yang, Wu, Li and Gao2019). Combined with the multiphase LB model and the enthalpy thermal LB model, Zhang et al. (Reference Zhang, Zhang, Fang, Zhao and Yang2020) proposed an axisymmetric LB model to simulate the freezing process of a sessile water droplet with consideration of droplet volume expansion. The combination of the pseudo-potential multiphase model and a thermal single-component phase change model is used by Pérez et al. (Reference Pérez, Leclaire, Ammar, Trépanier, Reggio and Benmeddour2021) for studying droplet impact and solidification on an airfoil.
Although multiphase flows have been investigated widely with different LB schemes, thermal multiphase simulations of liquid–solid phase transition and possible volume change during this transition (Gong et al. Reference Gong, Hou, Yang, Wu, Li and Gao2019; Zhang et al. Reference Zhang, Zhang, Fang, Zhao and Yang2020; Pérez et al. Reference Pérez, Leclaire, Ammar, Trépanier, Reggio and Benmeddour2021) are almost exclusively limited to the pseudo-potential approach. The method's limitation is that it relies on mesoscale interaction to separate multiple phases, making it difficult to extend to other (e.g. CFD) methods. On the other hand, the free energy approach uses novel mathematical physics to explain the phase separation process. This concept, referred to as the Cahn–Hilliard equation on a macro scale, can be easily understood on a mesoscopic level by utilizing the LB equation. Therefore, any improvements made to the free energy method can be effortlessly applied to the mesoscopic and macroscopic CFD methods.
As far as the authors know, there is no known free energy method for dealing with volume changes in multiphase flow. To fill this gap, this paper aims to present a new modification to the Cahn–Hilliard equation that can simulate volume changes during solidification. Unlike previous methods (Gong et al. Reference Gong, Hou, Yang, Wu, Li and Gao2019; Zhang et al. Reference Zhang, Zhang, Fang, Zhao and Yang2020; Pérez et al. Reference Pérez, Leclaire, Ammar, Trépanier, Reggio and Benmeddour2021), this modification can be used for both expansion and contraction scenarios, not just expansion. To implement this new formulation, a new LB framework has been developed to handle property variation that can occur during solidification. To show the capability of the introduced scheme, here we focus on a practical problem: the freezing of a sessile droplet.
The freezing of a liquid droplet is of interest for many engineering processes such as spray coating, laser welding, additive manufacturing, icing on aircraft surfaces and others. Consequently, the solidification of droplets has been studied extensively over the past few decades. Efforts have been made to understand the droplet freezing process through theoretical models (Schultz, Worster & Anderson Reference Schultz, Worster and Anderson2001; Schetnikov, Matiunin & Chernov Reference Schetnikov, Matiunin and Chernov2015; Zhang et al. Reference Zhang, Wu, Min and Liu2017b), experiments (Jung et al. Reference Jung, Tiwari, Doan and Poulikakos2012; Zhang, Wu & Min Reference Zhang, Wu and Min2017a; Shi & Duan Reference Shi and Duan2022) or numerical simulations (Vu Reference Vu2018; Zhang et al. Reference Zhang, Liu, Wu and Min2018; Bodaghkhani & Duan Reference Bodaghkhani and Duan2020). As a numeric method, LB has also been used in a number of studies (Xiong & Cheng Reference Xiong and Cheng2018b,Reference Xiong and Chenga; Gong et al. Reference Gong, Hou, Yang, Wu, Li and Gao2019; Zhang et al. Reference Zhang, Zhang, Fang, Zhao and Yang2020; Pérez et al. Reference Pérez, Leclaire, Ammar, Trépanier, Reggio and Benmeddour2021). Most of these studies are concerned with the solidification of water droplets surrounded by cold air. Thus, the focus of the investigations was primarily on the droplet thermal conditions or wetting properties of the surface. But effects of the environmental (e.g. air) conditions are often ignored (Schetnikov et al. Reference Schetnikov, Matiunin and Chernov2015; Zhang et al. Reference Zhang, Liu, Wu and Min2018) or simplified to a constant temperature boundary condition at the droplet interface (Feuillebois et al. Reference Feuillebois, Lasek, Creismeas, Pigeonneau and Szaniawski1995; Tabakova & Feuillebois Reference Tabakova and Feuillebois2004). While these hypotheses may be true for water droplets in the presence of air, they cannot be applied to other applications, such as spray coating or gas shielded welding, where other materials may be involved. Nevertheless, these hypotheses can be challenged even with water droplets (Schultz et al. Reference Schultz, Worster and Anderson2001; Jung et al. Reference Jung, Tiwari, Doan and Poulikakos2012).
This study attempts to identify the individual and interacting effects of environmental conditions on the solidification process, which have been rarely investigated in previous studies. Furthermore, we examine the effects of positive as well as negative volume expansion rates on the final shape of solidified droplets. While droplet shape deformation due to volume expansion has been previously reported, no numerical simulation of solidification with volume contraction has been conducted before (including non-LB studies).
This paper is organized as follows. Section 2 briefly describes the sessile droplet solidification through macroscopic governing equations. Section 3 introduces the new LB framework. Boundary conditions and a discussion on choosing LB parameters are also included in this section. Next, in § 4 the introduced scheme is validated against classic numerical and experimental results available in the literature. Section 5 is devoted to numerical results, in which freezing time and shape deformation are investigated in detail. Finally, a conclusion is presented in § 6 to sum up the findings.
2. Macroscopic governing equations
This study focuses on a binary system including the ‘
$A$’ and the ‘
$B$’ components. For generality of our discussion, it is assumed that the 
$A$ component is going to retain its initial phase (liquid or gas) throughout the simulation. In contrast, the 
$B$ component experiences phase transitions (liquid to solid or vice versa) depending on the temperature conditions. The primary geometry is a liquid droplet (component 
$B$) settled on a cold surface whose temperature is kept constant at a value below the freezing temperature of the liquid. The liquid is surrounded by the 
$A$ component, whose temperature is also affected by the wall temperature. Since the droplet is static during solidification, inertia force is not important in the freezing process. In this regard, the density ratio was kept constant at a moderate range, e.g. 
${{{\rho _{{B}}}} / {{\rho _{{A}}}}} = 7$. Other parameters to consider are the substrate, droplet and environmental properties.
The cold wall temperature (
$T_{{wall}}$) and the contact angle (
$\psi _{{wall}}$) represent the substrate condition. Our analysis considered three thermal properties of the droplet, including its initial temperature (
${T_{{{ini(B)}}}}$) and its thermal diffusivity in solid (
${\alpha _{{{B(S)}}}}$) and liquid (
${\alpha _{{{B(L)}}}}$) states. Gravity (
$g$) is also included in this study as an independent variable. For environmental conditions, there are two important factors: initial temperature (
${T_{{{ini(A)}}}}$) and thermal diffusivity (
$\alpha _A$). It should be pointed out that in this definition of the problem, the droplet can be but is not necessarily water, and the surrounding fluid not necessarily air. They can be any two immiscible fluids whose conditions are determined by the aforementioned properties. Therefore, the results are applicable for a wide range of applications, such as welding with shielding gas, plasma-sprayed coatings and three-dimensional (3-D) metal printing. Figure 1 presents a schematic of the droplet and all parameters mentioned above.
Geometry, boundary condition (BC) and simulation parameters used for sessile droplet solidification.
In the free energy approach, two components are described by a numerical parameter – the component index 
$\phi$. This index generally varies between two specific values, e.g. 
$+\phi _0$ and 
$-\phi _0$ representing two bulk components. The variation of the component index between these two limiting values occurs through a diffuse interface that occupies a limited number of grid points. Here, we have
\begin{equation} {\rm{Composition}} = \left\{ {\begin{array}{@{}ll} {{A}}, & {\phi ={-} {\phi_0}}, \\ {{\rm{interface}}}, & {- {\phi_0} < \phi <{+} {\phi_0}}, \\ {{B}}, & {\phi ={+} {\phi_0}}. \end{array}}\right. \end{equation}The flow field of multicomponent Newtonian fluids with constant properties are governed by the continuity and momentum (the Navier–Stokes) equations
$$\begin{gather} {\partial_t}\rho + {\partial_\alpha }\left( {\rho {u_\alpha }} \right) = 0, \end{gather}$$
$$\begin{gather}\rho \left( {{\partial_t}{u_\alpha } + {u_\beta }{\partial_\beta }{u_\alpha }} \right) ={-} {\partial_\alpha }p + \mu {\partial_\beta }\left( {{\partial_\beta }{u_\alpha } + {\partial _\alpha }{u_\beta }} \right) + {F_\alpha }, \end{gather}$$
where 
${\boldsymbol {F}}={\boldsymbol {F}}_s+{\boldsymbol {F}}_g$ is the total body force including gravitational body force 
${\boldsymbol {F}}_g$ and the interfacial force 
${\boldsymbol {F}}_s$. Note that in all equations through this paper, the Einstein summation convention is considered for spatial coordinates, e.g. 
$\alpha$, 
$\beta$. For a binary system far from a solid surface, the equilibrium properties can be described by a Landau free energy functional of the form (Zheng et al. Reference Zheng, Shu and Chew2006)
\begin{equation} \varPsi = \int {\left( {\psi \left( \phi \right) + k\left( {{\partial_\alpha }\phi {\partial _\alpha }\phi } \right) + \frac{{{\rho_0}\ln {\rho_0}}}{3}} \right){\rm d}V}. \end{equation} The first term shows the bulk free energy 
$\psi (\phi ) = \mathcal {A}{({\phi ^2} - \phi _0^2)^2}$; the second term is responsible for the interaction between two fluids at the interface and the third term is the ideal gas free energy. The average density is denoted 
${\rho _0} = {{( {{\rho _{{A}}} + {\rho _{{B}}}} )} / 2}$. The term ‘
$\mathcal {A}$’ in the definition of free bulk energy and the term ‘
$k$’ used in the interaction term are two numerical parameters that control the surface tension 
$\sigma$ and width of the interface 
$\chi$,
$$\begin{gather} \sigma = \tfrac{4}{3}\sqrt {2k\mathcal{A}} \phi_0^3, \end{gather}$$
$$\begin{gather}\chi = \frac{1}{{{\phi_0}}}\sqrt {\frac{{2k}}{\mathcal{A}}}. \end{gather}$$The chemical potential is defined as the derivative of functional equation (2.4) with respect to phase index,
\begin{equation} {\mu_\phi } = 4\mathcal{A}( {{\phi ^3} - \phi_0^2\phi } ) - k{\partial_\alpha }{\partial_\alpha }\phi. \end{equation} The 
$\phi _0$ value is a positive non-zero value that logically can be chosen freely. Far from the interface where the second term in (2.7) is negligible, the chemical potential is equal for both components, i.e. 
${{\mu _{\phi ({{A}})}}={\mu _{\phi ({{B}})}}=0}$. Inside the interface the non-zero value of the first term balances the second term leading to again a zero chemical potential 
${\mu _{\phi ({\rm {interface}})}} = 0$. Any imbalance of chemical potential should lead to mass transfer to obtain a uniform distribution. That is why the interfacial body force in the Navier–Stokes equations (2.3) is defined as
\begin{equation} {F_{s,\alpha }} ={-} \phi {\partial_\alpha }{\mu_\phi }. \end{equation}For the component index, we use a modified version of the convective Cahn–Hilliard equation,
\begin{equation} {\partial_t}\phi + {\partial_\alpha }\left( {\phi {u_\alpha }} \right) = M{\partial_\beta }{\partial_\beta }{\mu_\phi } + S, \end{equation}
where ‘
$M$’ is the mobility. In the classical Cahn–Hilliard equation (Cahn & Hilliard Reference Cahn and Hilliard1958), interfacial diffusion fluxes are assumed to be proportional to chemical potential gradients. The equation is later extended to incorporate convective effects (Badalassi, Ceniceros & Banerjee Reference Badalassi, Ceniceros and Banerjee2003). The term (
$S$) is a fictional source term that we introduced to control volume variations during solidification. More details about this source term are provided in § 3.
Binary systems consist of two components, ‘A’ and ‘B’, which are tracked by the component index. It is assumed that the A component will not experience any phase change during the simulation. Component B, however, is going to solidify so it can be solid or liquid depending on its temperature. The solid phase can be described with two dimensionless numbers: the liquid fraction that shows what fraction of component B is liquid 
$(\,{{f_{{L}}} = {{{V_{{{B}}({{L}})}}} / {{V_{{B}}}}}} )$, and the solid fraction that indicates the amount of solid in the total volume 
$(\,{{f_{{S}}} = {{{V_{{{B}}({{S}})}}} / V}})$. These two fractions are related through the component index,
\begin{equation} {f_{{S}}} = \left( {1 - {f_{{L}}}} \right)\frac{{\phi + {\phi_0}}}{{2{\phi_0}}}. \end{equation}In this study, it is assumed that the thermodynamic properties of the B component vary during the solidification. In this regard, we use ‘B(S)’ and ‘B(L)’ to distinguish the B component properties in solid and liquid phases, respectively.
The temperature distribution is governed by the energy equation
\begin{equation} {d_t}H = {\partial_\beta }({k{\partial_\beta }T}). \end{equation}
The enthalpy ‘
$H$’ can be calculated as
\begin{equation} H = {H_{{A}}} + \frac{{\phi + {\phi_0}}}{{2{\phi_0}}}\left( {{H_{{B}}} - {H_{{A}}}} \right) = {(\rho c)_{{{eff}}}}T + \frac{{\phi + {\phi_0}}}{{2{\phi_0}}}{\rho _{{{B(L)}}}}{f_{{L}}}{L_{{h}}}, \end{equation}
where 
$L_{{h}}$ is the latent heat and 
$(\rho c)_{{{eff}}}$ is the effective volumetric heat capacity defined as
\begin{equation} {(\rho c)_{{{eff}}}} = \frac{{{\phi_0} - \phi }}{{2{\phi_0}}}{(\rho c)_{{A}}} + \frac{{{\phi _0} + \phi }}{{2{\phi_0}}}{(\rho c)_{{{B(L)}}}}{f_{{L}}} + {(\rho c)_{{{B(S)}}}}{f_{{S}}}. \end{equation}In a similar manner, the effective conductivity can be defined as a linear combination of component conductivities,
\begin{equation} {k_{{{eff}}}} = \frac{{{\phi_0} - \phi }}{{2{\phi_0}}}{k_{{A}}} + \frac{{{\phi_0} + \phi }}{{2{\phi_0}}}{k_{{{B(L)}}}}{f_{{L}}} + {k_{{{B(S)}}}}{f_{{S}}}. \end{equation}Substituting (2.12), (2.13) and (2.14) into (2.11) leads to
\begin{align} {d_t}T - \frac{{{k_{{{eff}}}}}}{{{{(\rho c)}_{{{eff}}}}}}{\partial_\beta }{\partial_\beta }T + \frac{{\phi + {\phi_0}}}{{2{\phi_0}}}\frac{{{\rho_{{{B(L)}}}}{L_{{h}}}}}{{{{(\rho c)}_{{{eff}}}}}}{d_t}{f_{{L}}} = \frac{{{\partial_\beta }T{\partial_\beta }{k_{{{eff}}}}}}{{{{(\rho c)}_{{{eff}}}}}} - \frac{{{f_{{L}}}{\rho_{{{B(L)}}}}{L_{{h}}}}}{{2{\phi_0}{{(\rho c)}_{{{eff}}}}}}{d_t}\phi - T{d_t}{(\rho c)_{{{eff}}}}. \end{align}The right-hand side of (2.15) is zero in every single-phase region (the A component and the B component in liquid and solid states) but non-zero at the interfaces. When the interface thickness is small compared with the characteristic length and the grid resolution is adequate, these terms can be ignored relative to the left-hand side of the equation. Thus, (2.15) will be simplified to
\begin{equation} {\partial_t}T + {\partial_\alpha }\left( {T{u_\alpha }} \right) = \frac{{{k_{{{eff}}}}}}{{{{(\rho c)}_{{{eff}}}}}}{\partial_\beta }{\partial_\beta }T - \frac{{\phi + {\phi_0}}}{{2{\phi_0}}}\frac{{{\rho_{{{B(L)}}}}{L_{{h}}}}}{{{{(\rho c)}_{{{eff}}}}}}{\partial_t}{f_{{L}}}. \end{equation}All governing equations ((2.2), (2.3), (2.9) and (2.16)) are solved using different LB schemes described in the next section.
3. Lattice Boltzmann method
In the LB method the physics of the problem is tracked through the evolution of fictitious particles named distribution functions, e.g. 
$m_i (x,t)$. At time 
$t$ for a grid point located at 
$x$, a series of distributions are defined corresponding to a set of discretized mesoscopic velocities indicated as a subscript index ‘
$i$’ (e.g. 
$m_1$, 
$m_2$). The lattice evolution equation for distribution functions is given below,
\begin{equation} {m_i}(x + {e_i}\delta t,t + \delta t) = {m_i}(x,t) + {\varOmega_{i,{{m}}}} + {S_{i,{{m}}}}, \end{equation}
where 
${\varOmega _{i,{{m}}}}$ is the collision term and 
${S_{i,{{m}}}}$ represents the source term or external body force effect if there was any. In the single relaxation time model the collision term is defined as follows:
\begin{align} {\varOmega_{i,{{m}}}} ={-} \frac{1}{{{\tau_{{m}}}}}\left( {{m_i}(x,t) - m_i^{eq}(x,t)} \right)\!. \end{align}
In this equation, 
$\tau _{{m}}$ is the dimensionless relaxation time and 
$m_i^{eq}$ is the equilibrium distribution. Generally, equilibrium distribution is defined as a function of macroscopic properties, mesoscopic velocities and a weighting factor.
There is a multi-scale analysis named Chapman–Enskog analysis that shows which macroscopic equation can be solved using the LB evolution equation. Note that the LB equation introduced in this paper is not just a flow solver, but a numerical tool to solve all governing equations ((2.2), (2.3), (2.9) and (2.16)). The solution of each governing equation requires at least three definitions: distribution functions, equilibrium distributions and dimensionless relaxation time.
In the following subsections, three sets of definitions are introduced to reproduce the governing equations mentioned in § 2. The LB equations introduced in this paper can be used (without any modifications) for two-dimensional (2-D) and 3-D simulations using D2Q9, D3Q19 or D3Q27 lattices. Table 1 provides a list of the mesoscopic velocities and weighting factors for these lattices.
Mesoscopic velocities and weighting factors for 2-D and 3-D lattices.
3.1. Multiphase flows
For simulating a binary system here, we used a modified version of what was proposed by Shao et al. (Reference Shao, Shu, Huang and Chew2014). According to this model, which is a recasting of He et al.'s approach (1999), the evolution equation for the density distribution is given by
\begin{equation} {f_i}(x + {e_i}\delta t,t + \delta t) = {f_i}(x,t) - \frac{1}{{{\tau_{{f}}}}}\left(\,{{f_i}(x,t) - f_i^{{{eq}}}(x,t)} \right) + {G_i} + {F_i}, \end{equation}
where 
$G_i$ represents the body force effects, including multiphase force 
$(-\phi \partial _\alpha \mu _\phi )$ and gravitational body forces 
$(\rho g_\alpha )$, and 
$F_i$ shows the pressure effects due to different densities (Shao et al. Reference Shao, Shu, Huang and Chew2014). These two terms are defined as follows:
$$\begin{gather} {G_i} = \left( {1 - \frac{1}{{2{\tau_{\rm{f}}}}}} \right)\delta t{w_i}\left( {1 + \frac{{{e_{i\beta }}{u_\beta }}}{{c_{{s}}^2}} + \frac{{{e_{i\beta }}{e_{i\gamma }}{u_\beta }{u_\gamma }}}{{2c_{{s}}^4}} - \frac{{{u_\beta }{u_\beta }}}{{2c_{{s}}^2}}} \right)\left( {\frac{{{e_{i\alpha }} - {u_\alpha }}}{{c_{{s}}^2}}} \right)( {\rho {g_\alpha } - \phi {\partial_\alpha }{\mu_\phi }}), \end{gather}$$
$$\begin{gather}{F_i} = \left( {1 - \frac{1}{{2{\tau_{{f}}}}}} \right)\delta t{w_i}\left( {\frac{{{e_{i\beta }}{u_\beta }}}{{c_{{s}}^2}} + \frac{{{e_{i\beta }}{e_{i\gamma }}{u_\beta }{u_\gamma }}}{{2c_{{s}}^4}} - \frac{{{u_\beta }{u_\beta }}}{{2c_{{s}}^2}}} \right)\left( {\frac{{{e_{i\alpha }} - {u_\alpha }}}{{c_{{s}}^2}}} \right)\left( {c_{{s}}^2{\partial_\alpha }\rho } \right)\!. \end{gather}$$The force term formulation in (3.4) and (3.5) is based on a discrete Boltzmann equation proposed by He, Shan & Doolen (Reference He, Shan and Doolen1998b). Guo, Zheng & Shi (Reference Guo, Zheng and Shi2002) explained how this formulation may not be consistent with the hydrodynamic equations, particularly where the body force changes both in time and space. Inspired by the idea used in Guo et al. (Reference Guo, Zheng and Shi2002), we redesigned the forcing terms to apply body force and pressure difference effects more accurately and efficiently, i.e.
$$\begin{gather} G_i^{{{new}}} = \left( {1 - \frac{1}{{2{\tau_{{f}}}}}} \right)\delta t{w_i}\left( {\frac{{{e_{i\alpha }} - {u_\alpha }}}{{c_{{s}}^2}} + \frac{{{e_{i\beta }}{u_\beta }{e_{i\alpha }}}}{{c_{{s}}^4}}} \right)( {\rho {g_\alpha } - \phi {\partial_\alpha }{\mu_\phi }} ), \end{gather}$$
$$\begin{gather}F_i^{{{new}}} = \left( {1 - \frac{1}{{2{\tau_{{f}}}}}} \right)\delta t{w_i}\left( {\frac{{{e_{i\beta }}{u_\beta }{e_{i\alpha }}}}{{c_{{s}}^4}}} \right)\left( {c_{{s}}^2{\partial_\alpha }\rho } \right)\!. \end{gather}$$Equilibrium distribution in (3.3) is expressed as
\begin{equation} f_i^{{{eq}}}(x,t) = {w_i}\rho \left( {\frac{{{\rho_0}}}{\rho } + \frac{{{e_{i\beta }}{u_\beta }}}{{c_{{s}}^2}} + \frac{{{e_{i\beta }}{e_{i\gamma }}{u_\beta }{u_\gamma }}}{{2c_{{s}}^4}} - \frac{{{u_\beta }{u_\beta }}}{{2c_{{s}}^2}}} \right)\!. \end{equation}The macroscopic properties can be calculated as
$$\begin{gather} \rho = {\rho_{{A}}} + \frac{{\phi + {\phi_0}}}{{2{\phi_0}}}\left( {{\rho_{{B}}} - {\rho_{{A}}}} \right)\!, \end{gather}$$
$$\begin{gather}{\rho_0} = \sum_i {{f_i}} + \frac{1}{2}\delta t\left( {{u_\alpha }{\partial_\alpha }\rho } \right)\!, \end{gather}$$
$$\begin{gather}\rho {u_\alpha } = \sum_i {{f_i}{e_{i\alpha }}} + \frac{1}{2}\delta t( {\rho {g_\alpha } - \phi {\partial_\alpha }{\mu_\phi }} ). \end{gather}$$ The dimensionless relaxation time is related to the dynamic viscosity through 
$\mu = \rho c_{{s}}^2\delta t( {{\tau _{{f}}} - 0.5} )$. Far from the interface, this equation can be used for each fluid where each of them has its own dimensionless relaxation time, i.e.
$$\begin{gather} {\mu_{{A}}} = {\rho_{{A}}}c_{{s}}^2\delta t( {{\tau_{{{f(A)}}}} - 0.5}), \end{gather}$$
$$\begin{gather}{\mu_{{B}}} = {\rho_{{B}}}c_{{s}}^2\delta t( {{\tau_{{{f(B)}}}} - 0.5} ). \end{gather}$$ For the interface, it is common to use a linear approximation based on the component index value 
$\mu ( \phi ) = {\mu _{{A}}} + ( {{\mu _{{B}}} - {\mu _{{A}}}} ){{( {\phi + {\phi _0}} )} / {2{\phi _0}}}$ (Kuzmin et al. Reference Kuzmin, Januszewski, Eskin, Mostowfi and Derksen2011). Using this equation, the dimensionless relaxation time can be obtained as
\begin{equation} {\tau_{{f}}} = \frac{{\mu \left( \phi \right)}}{{\rho c_{{s}}^2\delta t}} + \frac{1}{2}, \end{equation}
where 
$\rho$ is the local density obtained from (3.9). Like all other schemes designed using He et al.'s approach (1999), the model is restricted to moderate density ratios. In static cases or situations with negligible inertia effects, a moderate density ratio or even a density-matched model can be safely used (Kuzmin et al. Reference Kuzmin, Januszewski, Eskin, Mostowfi and Derksen2011). When simulating multiphase flows with high-density ratios is inevitable, an alternative, more computationally expensive model (Inamuro et al. Reference Inamuro, Ogata, Tajima and Konishi2004; Lee & Lin Reference Lee and Lin2005) can be used to replace (3.3).
The dynamics of the fluid interface is governed by the modified convective Cahn–Hilliard equation (2.9). This equation can be solved numerically using the following LB equation:
\begin{equation} {h_i}(x + {e_i}\delta t,t + \delta t) = {h_i}(x,t) - \frac{1}{{{\tau_{{h}}}}}\left( {{h_i}(x,t) - h_i^{{{eq}}}(x,t)} \right) + {S_{i,{h}}}. \end{equation}The equilibrium distribution and the dimensionless relaxation time are defined as
$$\begin{gather} h_i^{{{eq}}}(x,t) = \left\{ \begin{array}{@{}ll} {{w_i}\left( {\dfrac{{\varGamma {\mu_\phi }}}{{c_{{s}}^2}} + \phi \dfrac{{{e_{i\beta }}{u_\beta }}}{{c_{{s}}^2}} + \phi \dfrac{{{e_{i\beta }}{e_{i\gamma }}{u_\beta }{u_\gamma }}}{{2c_{{s}}^4}} - \phi \dfrac{{{u_\beta }{u_\beta }}}{{2c_{{s}}^2}}} \right)}, & {{{i}} > 0}, \\[12pt] {\phi - \sum\nolimits_{i \ne 0} {h_i^{{{eq}}}} }, & {{{i}} = 0,} \end{array} \right. \end{gather}$$
$$\begin{gather}{\tau_{{h}}} = \frac{M}{{\varGamma \delta t}} + \frac{1}{2}, \end{gather}$$
where 
$\varGamma$ is a numerical parameter. The component index value can be calculated through
\begin{equation} \phi = \sum_i {{h_i}}. \end{equation} In the LB scheme a source term 
$S_{i,{h}}$ can be implemented using the following equation (Shi & Guo Reference Shi and Guo2009):
\begin{equation} {S_{i,{h}}} = {w_i}S\left( {1 + \left( {1 - \frac{1}{{2{\tau_{{h}}}}}} \right)\frac{{{e_{i\alpha }}{u_\alpha }}}{{c_{{s}}^2}}} \right)\!. \end{equation} The source value ‘
$S$’ in the above equation is the new source term we introduced in (2.9) to control volume changes during solidification. As mentioned earlier, components are tracked using the component index. As a result, any change in the volume of a component must be represented by a change in its index value. Here the idea is to manipulate the component index value to simulate the volume change of a component. By modifying the component index value, the system deviates from its equilibrium condition, resulting in a non-zero chemical potential value. Consequently, the diffusion mechanism will be activated, balancing the chemical distribution through expansion and compression (corresponding to positive and negative volume expansion values). The density values are then modified in the energy equation (2.16) to conserve the total mass. The key question is to what extent the value of the index should change for a given volume expansion rate.
The share of one component (here the B component) from a typical volume of the computational domain can be described as 
${V_{{B}}} = V{{(\phi + {\phi _0})} / {(2{\phi _0})}}$. By differentiating this equation, it is possible to find the 
$\phi$ modification required for a given volume expansion, i.e.
\begin{equation} \delta {V_{{B}}} = \frac{{\delta \phi }}{{2{\phi_0}}}V = \frac{{\delta \phi }}{{2{\phi_0}}}\left( {\frac{{2{\phi_0}}}{{\phi + {\phi_0}}}{V_B}} \right) = \frac{{\delta \phi }}{{\phi + {\phi_0}}}{V_{{B}}}. \end{equation}During the solidification, the volume of the B component will change, i.e.
\begin{equation} \delta {V_{{B}}} = \left( {{V_{{{B(L)}}}} + {f_{{S}}}\left( {{V_{{{B(S)}}}} - {V_{{{B(L)}}}}} \right)} \right) - {V_{{{B(L)}}}}. \end{equation} Using (3.20) and (3.21), it can be concluded that, for a given expansion rate such as 
$\lambda = {{({V_{{{B(S)}}}} - {V_{{{B(L)}}}})} / {{V_{{{B(L)}}}}}}$, the necessary change of the index value would be
\begin{equation} \delta \phi = \lambda \left( {\phi + {\phi_0}} \right)\delta {f_{{S}}}. \end{equation} It is important to note that solidification is not an instantaneous process, and (3.22) is going to apply step by step during the time-marching solution of the governing equations. This means that for a constant 
$\lambda$ value, the final volume will be obtained using a geometric progression
\begin{equation} {V_{{{final}}}} = \sum_{k = 0}^\infty {{V_{{{initial}}}}{\lambda ^k}} = {V_{{{initial}}}}\frac{{1 - {\lambda ^\infty }}}{{1 - \lambda }} = \frac{{{V_{{{initial}}}}}}{{1 - \lambda }}. \end{equation}Considering this point in the calculations, we introduce the new sourcing term for the Cahn–Hilliard equation as
\begin{equation} S = \left( {\phi + {\phi_0}} \right)\frac{\lambda }{{1 + \lambda }}{\partial_t}{f_{{S}}}. \end{equation} By this definition, a complete solidification theoretically should lead to 
${V_{{{initial}}}}( {1 + \lambda } )$ as the final volume. But in practice, the complexity of multiphase flows may deviate the final volume from what is expected. More discussion about this deviation is provided in § 3.5.
3.2. Thermal simulation
The LB evolution equation for the temperature field is given by
\begin{equation} {l_i}(x + {e_i}\delta t,t + \delta t) = {l_i}(x,t) - \frac{1}{{{\tau_{{l}}}}}\left( {{l_i}(x,t) - l_i^{{{eq}}}(x,t)} \right) + {S_{i,{{l}}}}, \end{equation}where the equilibrium distribution and the source term are defined as
$$\begin{gather} l_i^{{{eq}}}(x,t) = {w_i}T\left( {1 + \frac{{{e_{i\beta }}{u_\beta }}}{{c_{{s}}^2}} + \frac{{{e_{i\beta }}{e_{i\gamma }}{u_\beta }{u_\gamma }}}{{2c_{{s}}^4}} - \frac{{{u_\beta }{u_\beta }}}{{2c_{{s}}^2}}} \right)\!, \end{gather}$$
$$\begin{gather}{S_{i,{{l}}}} = -{w_i}\left( {1 + \left( {1 - \frac{1}{{2{\tau_{{l}}}}}} \right)\frac{{{e_{i\alpha }}{u_\alpha }}}{{c_{{s}}^2}}} \right)\left( {\frac{{\phi + {\phi_0}}}{{2{\phi_0}}}\frac{{{\rho_{{{B(L)}}}}{L_{{h}}}}}{{{{\left( {\rho c} \right)}_{{{eff}}}}}}{\partial_t}{f_{{L}}}} \right)\!. \end{gather}$$The temperature value at each grid location is obtained as
\begin{equation} T = \sum_i {{l_i}}. \end{equation}The dimensionless relaxation time is defined based on local effective properties provided in (2.13), (2.14),
\begin{equation} {\tau_{{l}}} = \frac{{{k_{{{eff}}}}}}{{{{(\rho c)}_{{{eff}}}}}}\frac{1}{{c_{{s}}^2\delta t}} + \frac{1}{2}. \end{equation}3.3. Solid domain consideration
During the simulations some part of the computational domain may solidify. For these parts, the distributions must remain in the solid domain. This can be applied by changing the lattice evolution equation (3.1) into
\begin{equation} {m_i}(x + {e_i}\delta t,t + \delta t) = {m_i}(x,t) + ( {1 - {B_{{m}}}} )( {{\varOmega_{i,{{m}}}} + {S_{i,{{m}}}}} ) + {B_{{m}}}\varOmega_{i,{{m}}}^{{{solid}}}. \end{equation} In this equation, 
$B_{{m}}$ is a numerical parameter whose value varies between 
$0$ to 
$1$ depending on the fluid fraction value. The new term 
$\varOmega _{i,{{m}}}^{{{solid}}}$ is responsible to keep information in the solid domain. For the density distribution functions, the approach provided in Huang et al. (Reference Huang, Wu and Cheng2013) is used as follows:
$$\begin{gather} {B_{{f}}} = \frac{{{\tau_{{f}}} - 0.5}}{{{\tau_{{f}}} - 0.5 + {f_{{L}}}}}\left( {1 - {f_{{L}}}} \right)\!, \end{gather}$$
$$\begin{gather}\varOmega_{i,{{f}}}^{{{solid}}} = f_i^{{{eq}}({{solid}})} + {f_{{{opp}}(i)}}(x,t) - f_{{{opp}}(i)}^{{{eq}}} - {f_i}(x,t). \end{gather}$$ In this equation, 
${\rm {opp}}(i)$ shows the direction opposite to the 
$i$ direction. The term 
$f_i^{{{eq}}({{solid}})}$ is the solid equilibrium distribution and can be calculated using (3.8), but in this calculation the velocity should be set equal to the solid phase velocity that in this case is zero.
For component index distributions, we introduce the following definitions:
$$\begin{gather} {B_{{h}}} = 1 - {f_{{L}}}, \end{gather}$$
$$\begin{gather}\varOmega_{i,{{h}}}^{{{solid}}} = {h_{{{opp}}(i)}}(x,t) - {h_i}(x,t). \end{gather}$$ Note that the volume change during solidification depends on the partial differential of the liquid fraction with respect to time. When the liquid fraction is lower than unity and going to increase in the following simulation steps, the solid collision may trap the produced phase index value inside the solid domain. To avoid this situation, the solid collision (3.34) must be applied only for nodes with constant liquid fractions 
$( {{\partial _t}{f_{{L}}} = 0} )$.
The temperature LB evolution equation does not need to use a solid collision function because the no-slip velocity in the solid domain is actually applied using (3.32).
3.4. Boundary conditions
Three boundary condition schemes are used to handle wall physics. In the flow field one of the most common methods to apply a no-slip boundary condition is a so-called ‘bounce-back rule’. Despite its simplicity and ease of implementation, it may lead to an undesirable slip velocity depending on the viscosity value (He et al. Reference He, Zou, Luo and Dembo1997; Mohammadipoor, Niazmand & Mirbozorgi Reference Mohammadipoor, Niazmand and Mirbozorgi2014). The problem is addressed before by introducing a general boundary treatment (Mohammadipour, Niazmand & Succi Reference Mohammadipour, Niazmand and Succi2017). The wettability of the solid surface is modelled using (Briant & Yeomans Reference Briant and Yeomans2004). The boundary condition in the Cahn–Hilliard equation is not limited to the wettability condition. The preservation of mass in the computational domain requires us to apply zero normal gradients for chemical potential 
$( {{{{{\partial _n}\mu } |}_{{{wall}}}} = 0} )$ at the boundaries (Li et al. Reference Li, Luo, Kang, He, Chen and Liu2016). This equation must be considered in the collision step for boundary nodes where the interfacial force is calculated using (2.8).
For the energy equation, the desired boundary temperature can be applied to the boundary nodes using the following equation (Mohammadipour, Succi & Niazmand Reference Mohammadipour, Succi and Niazmand2018):
\begin{align} {l_{i({{unknown}})}}({x_{{{wall}}}},t) = \frac{{{w_i}}}{{\sum\nolimits_{{{unknown}}} {{w_i}} }}\left( {{T_{{{wall}}}} - \sum\limits_{{{known}}} {{l_i}({x_{{{wall}}}},t)} } \right)\!. \end{align}3.5. Selecting LB parameters
According to the similarity concept, the LB results have similitude with the real application as long as the two models share similar geometries and equal effective dimensionless numbers. For multiphase solidification simulated in this paper, the important dimensionless parameters are the Stefan number, the Atwood number and the Bond number, defined respectively as
$$\begin{gather} {{Ste}} = \frac{{c\left( {{T_{{m}}} - {T_{{{wall}}}}} \right)}}{{{L_h}}}, \end{gather}$$
$$\begin{gather}{{Atw}} = \frac{{{\rho_{{B}}} - {\rho_{{A}}}}}{{{\rho_{{B}}} + {\rho_{{A}}}}}, \end{gather}$$
$$\begin{gather}{{Bo}} = \frac{{\left( {{\rho_{{B}}} - {\rho_{{A}}}} \right)gl_{{{ref}}}^2}}{\sigma }. \end{gather}$$ There are also some numerical parameters such as 
$\phi _0$, 
$\varGamma$, 
$\mathcal {A}$,
$\ldots$ that do not have a clear counterpart of thermodynamics properties. While there is no solid rule for how these parameters should be set, it is important to note that they are strongly correlated with simulation accuracy and stability. Here we propose a step-by-step process for establishing simulation parameters.
First, the density distribution is initialized using the Atwood number and the volume expansion coefficient 
$\lambda$,
\begin{align} \rho = \left\{ {\begin{array}{{@{}ll}} {1 + {Atw}}, & {{B\ {\rm component\ in\ liquid\ phase}}}, \\ {{{(1 + {Atw})} / {(1 + \lambda )}}}, & {{B\ {\rm component\ in\ solid\ phase}}}, \\ {1 - {Atw}}, & {{A\ {\rm component}}}. \end{array}} \right. \end{align} It is very common to define the 
$\phi _0$ parameter as 
${\phi _0} = {{({\rho _{{B}}} - {\rho _{{A}}})} / 2}$ (Shao et al. Reference Shao, Shu, Huang and Chew2014), but since 
$\phi _0$ cannot be zero, the usability of this definition is limited to components with different densities. As a remedy, here we define 
$\phi _0$ in a new way, i.e.
\begin{equation} {\phi_0} = \frac{{2{\rho_{{A}}}}}{{{\rho_{{A}}} + {\rho_{{B}}}}}. \end{equation} The surface tension is considered as an input parameter and the interface thickness is set to 
$\chi =4 \sim 5$ as is suggested in Zheng et al. (Reference Zheng, Shu and Chew2006). Here ‘
$\mathcal {A}$’ and ‘
$k$’ parameters in the chemical potential definition are set as
$$\begin{gather} \mathcal{A} = \frac{3}{4}\frac{\sigma }{{\phi_0^4\chi }}, \end{gather}$$
$$\begin{gather}k = \tfrac{1}{2}\mathcal{A}{\chi ^2}\phi_0^2. \end{gather}$$ The dimensionless relaxation time for the Cahn–Hilliard equation is set to be 
$\tau _{{h}}=1$ to minimize unphysical spurious velocities that may appear close to interfaces (Pooley, Kusumaatmaja & Yeomans Reference Pooley, Kusumaatmaja and Yeomans2008), whereas other dimensionless parameters are chosen by thermodynamic properties of the A and B components. The heat capacity 
$(c)$, the latent heat 
$(L_h)$, the melting temperature 
$(T_{{m}})$ and the wall temperature 
$(T_{{wall}})$ are determined by the Stefan number definition (3.36).
For a chosen 
$\tau _{{h}}$, the 
$\varGamma$ value controls the mobility factor ‘
$M$’ in the Cahn–Hilliard equation. It should be set high enough to ensure that the interface has enough time to relax to local equilibrium but not too high that residual diffusion becomes significant (Kendon et al. Reference Kendon, Cates, Pagonabarraga, Desplat and Bladon2001). Here however there is another concern. The partial time derivative of the solid factor 
$f_{{S}}$ used in the source definition (3.24) is determined from the energy equation. This relation couples the energy equation to the Cahn–Hilliard equation and makes the volume change simulation a transient approach. For LB simulations where the time step for all equations is set to be one in the lattice unit system, the reference or the characteristic time becomes important. In the energy equation, a simple scale analysis between the diffusion and source terms shows that the time scale is
\begin{equation} {{{t_{{{ref}}}}} |_T} = \frac{{l_{{{ref}}}^2}}{{{\alpha_T}{{Ste}}}}, \end{equation}
where 
$\alpha _T$ is thermal diffusivity. The characteristic time reference for the Cahn–Hilliard equation can be obtained by balancing the transient and diffusion terms,
\begin{equation} { {{t_{{{ref}}}}} |_\phi } = \frac{{l_{{{ref}}}^3\phi_0^2}}{{\sigma {M}}}. \end{equation}A new dimensionless number can be defined as the ratio of these two time scales,
\begin{equation} {R_t} = \frac{{{{ {{t_{{{ref}}}}} |}_\phi }}}{{{{ {{t_{{{ref}}}}} |}_T}}} = \frac{{{l_{{{ref}}}}\phi_0^2{\alpha_T}{{Ste}}}}{{M\sigma }}. \end{equation} By choosing the interface thickness as the reference length, 
$R_t$ will be used to determine the mobility factor
\begin{equation} M = \frac{{\chi \phi_0^2{\alpha_T}{{Ste}}}}{{{R_t}\sigma }}. \end{equation} Numerical simulation shows that a balance between accuracy and computational cost can be achieved by keeping 
$R_t$ close to unity.
4. Validations
In this section the accuracy and stability of the developed framework are evaluated by simulating some classical multiphase problems including the Laplace test, the Rayleigh–Taylor instability, expansion and contraction of a droplet when it solidifies inside a stationary flow and freezing of a sessile water droplet on a cold surface.
4.1. The Laplace test
The Laplace test is about the equilibrium condition of a liquid drop or a gas bubble when it is surrounded by another fluid. In this condition the pressure difference between the inside and outside of the drop or the bubble is related to the surface tension and curvature of the interface. In a 2-D situation where the radius of the drop is 
$R$, the Laplace equation reads as 
$\Delta p = {\sigma / R}$, where 
$\sigma$ is the surface tension.
To simulate this test case, a square domain with the size of 
$6R \times 6R$ is considered where a circular drop is located at the centre of the computational domain. There is no gravity, all boundaries are periodic and the drop is modelled using 
$R=20$, 
$30$ and 
$40$ lattice grids. In the free energy approach the pressure tensor is obtained as (Kusumaatmaja & Yeomans Reference Kusumaatmaja and Yeomans2010)
\begin{equation} {P_{\alpha \beta }} = \left( {\mathcal{A}\left[ {3{\phi ^4} - 2\phi_0^2{\phi ^2} - \phi_0^4} \right] - k\phi {\partial_\gamma }{\partial_\gamma }\phi - \frac{k}{2}{\partial_\gamma }\phi {\partial_\gamma }\phi + \frac{{{\rho_0}}}{3}} \right){\delta_{\alpha \beta }} + k{\partial_\alpha }{\partial_\beta }\phi. \end{equation} The static pressure difference used in the Laplace equation is determined by calculating the scalar portion of the pressure tensor in the centre of the drop (for internal pressure) and at domain boundaries far from the drop (for external pressure). Figure 2 shows the variation of pressure difference as a function of the curvature 
$({1 / R})$. As is shown for increased grid resolution (reduced curvature value), the pressure difference obtained from the numerical simulations coincides with what is predicted by the Laplace equation.
Variation of pressure difference as a function of the curvature 
$(1/R)$ for two different values of surface tension.
4.2. The Rayleigh–Taylor instability
The Rayleigh–Taylor instability problem tracks the interface of two immiscible liquids with different densities where the lighter one is placed beneath the heavier one. A small perturbation at the interface will intensify by gravity leading to rising bubbles of lighter liquid and spikes of the heavier one falling. The flow is characterized by two dimensionless numbers: the Atwood number and the Reynolds number. The former is introduced before in § 3, and the latter is defined as
\begin{equation} {{Re}} = \frac{{\sqrt {Lg} L}}{\upsilon }, \end{equation}
where 
$L$ denotes the domain's width, 
$g$ is the gravity and 
${\sqrt {Lg} }$ represents the characteristic velocity. To simulate this test case, a rectangular domain with the size of 
$L \times 4L$ is used. The Atwood number is set to be 
${Atw}=0.5$, while the gravity is chosen in such a way that the Reynolds number is 
$256$. The simulation was conducted using two grid resolutions of 
$L=100$ and 
$200$. Initially, half of the domain is filled with light liquid with a density of 
$\rho _A=1-Atw=0.5$ and the upper half is set to have 
$\rho _B=1+{Atw}=1.5$. To simulate a perturbation, the interface is set to follow the equation
\begin{equation} y = 0.1L\cos \left( {\frac{{2{\rm \pi} x}}{L}} \right)\!. \end{equation}
Figure 3(a) shows the vertical position of the interface for the bubble and the spike during the time evolution of simulation and compares it with what is reported in He et al. (Reference He, Chen and Zhang1999). The time in this simulation is normalized as 
${{t^*}= {t / {\sqrt {{L / g}} }}}$. The flow configurations in three dimensionless times of 
$t^*=1$, 
$3$ and 
$5$ are depicted in figure 3(b). The results obtained in this study appears to compare well with those in He et al. (Reference He, Chen and Zhang1999).
(a) The vertical position of the interface for the bubble and the spike versus time: results of current study (open diamonds, and filled circles) and those of He et al. (Reference He, Chen and Zhang1999) (the solid line). (b) The evolution of fluid interface in the Rayleigh–Taylor instability at 
${Re} = 256$ where the grid resolution is 
$L=200$.
4.3. Volume change during solidification
This subsection evaluates how the new scheme can simulate the volume change when the volume expansion coefficient is positive or negative. The test case designed for this purpose is a drop of the B component surrounded by the A component inside a square domain (figure 4a). The drop and surrounding fluid are initially at a temperature higher than the melting temperature. The lower boundary is a cold wall with a temperature lower than the melting temperature, and its value is set by the Stefan number. The left and right boundaries are periodic, and a zero gradient of temperature is applied to the upper boundary. There is no gravity, so the drop remains stationary at its initial location. The drop volume is calculated before and after solidification using the following equation:
\begin{equation} V = \sum_{{\rm{all}}\ {\rm{nodes}}} {\frac{{\phi (x,y) + {\phi_0}}}{{2{\phi_0}}}}. \end{equation}(a) Simulation configuration for testing the volume change during solidification. (b) Volume change obtained from simulations versus nominal expansion rate values for three sets of simulations: 
${Atw}=0.75$, 
$D=40$, 
${Ste}=0.1$ (
$\vartriangle$), 
${Atw}=0.5$, 
$D=50$, 
${Ste}=0.125$ (
$\square$), and 
${Atw}=0.0$, 
$D=60$, 
${Ste}=0.15$ (
$\triangledown$).
For three density ratios (
${Atw}=0$, 
$0.5$ and 
$0.75$), three droplet sizes (
$D=40$, 
$50$ and 
$60$) and three Stefan numbers (
${Ste}=0.1$, 
$0.125$ and 
$0.15$), the change in volume obtained from the simulation (
$\lambda _{{numerical}}$) is shown in figure 4(b) with respect to the nominal expansion rate value (
$\lambda$) used in (3.24). The numerical results indicate that the introduced method can successfully simulate the volume change during a solidification. Figure 4(b) shows some deviation from expectations when the magnitude of volume change is high (
$| \lambda | > 10\,\%$). This is related to the complexity and transient behaviour of the volume change that was addressed before in § 3.5. Note that this deviation can be corrected by adjusting the input 
$\lambda$ value in (3.24) using the shooting method.
4.4. Water droplet freezing
Like any other substance, water has a specific freezing temperature for a given environmental pressure. Below this melting temperature, the stable form of water is solid. However, experimental results indicate that water can experience liquid form even at a temperature range lower than the melting temperature. This supercooled water is metastable and needs a minor perturbation or temperature fluctuation to initiate the nucleation process. The next step is the recalescence stage when heterogeneous or homogeneous nucleation converts the liquid water into a uniform mixture of ice and water. The temperature in this stage rises back to freezing temperature. Since the recalescence stage is very fast (it takes only a few milliseconds), the solidification in this stage is limited to a small portion of the droplet. After recalescence, the water–ice mixture, which has almost a uniform temperature, starts to solidify. The solidification process occurs at a constant temperature (melting temperature) and is accompanied by changes in water properties. Table 2 lists the thermal properties of ice and water.
Thermal properties of ice and water.
Here, we ignore the recalescence stage and use it only to initialize our numerical set-up. The simulation utilizes a 3-D domain with 19 discrete velocities at each grid point (D3Q19 lattice). The B component acts as a water drop, and its thermal properties are obtained from the water–ice mixture. The liquid mass fraction of water in the mixture can be approximated as (Blake et al. Reference Blake, Thompson, Raps and Strobl2014)
\begin{equation} \xi = 1 - \frac{{{{c}_{{{water}}}}}}{{{{c}_{{{ice}}}}}}{{Ste}}. \end{equation} Having the water mass fraction, any thermal property, e.g. 
$M_{{mix}}$ of the mixture, can be determined, i.e.
\begin{equation} {M_{{{mix}}}} = {M_{{{water}}}}\xi + {M_{{{ice}}}}\left( {1 - \xi } \right)\!. \end{equation} The experimental results used as the validation reference come from Shi & Duan (Reference Shi and Duan2022), where they examined the anti-icing performance of stainless-steel-based surfaces. The data we received from the authors is related to the freezing process of a water droplet with a volume of 
$14.76\ {\rm {mm}}^3$ on a cold surface with a static contact angle of 
$\psi _{{wall}}=106.86^\circ$. The wall temperature in this experiment was kept constant at 
$-14.6\,^\circ {\rm {C}}$. The dimensionless numbers used to conduct the simulation are listed in table 3.
Essential dimensionless parameters associated with water droplet freezing.
$^{a}$ The characteristic length in the Bond number definition is chosen to be 
$l_0=V/A$.
According to the experimental data, the solidification process took 
$t_{{freezing}} = 29.68$ s. Using this value as a reference time, the dimensionless time can be defined as 
$t^*=t/t_{{freezing}}$. Figure 5 compares the numerical results obtained from the LB simulation with those from the experiment at four different times 
$t^*=0$, 
$0.25$, 
$0.5$ and 
$1$ (images printed with authors’ permission). As is shown in figure 5, the numerical results compare well with the experimental data in contact angle (figure 5a), the shape and the evolution of the water–ice interface during the solidification (figure 5b,c). When water solidifies, its density decreases, and the increase in volume may lead to a conical pointy tip at the top of the ice. Figure 5(d) shows that the current LB scheme can successfully capture this phenomenon. The tip angle obtained from the simulation (
$127^\circ$) is in good agreement with the experimental data and what is reported in the literature (Marín et al. Reference Marín, Enríquez, Brunet, Colinet and Snoeijer2014; Schetnikov et al. Reference Schetnikov, Matiunin and Chernov2015).
Experimental and numerical droplet shape at (a) 
$t^*=0$, (b) 
$t^*=0.25$, (c) 
$t^*=0.5$ and (d) 
$t^*=1.0$.
5. Results and discussions
To make a droplet solidify, thermal energy must be extracted from the liquid. The heat transfer in this study is limited to conduction and the cold reservoir is simulated as a substrate with a constant temperature below the melting point of the liquid. The rate of heat transfer depends on the substrate and droplet properties introduced in § 2. The freezing process of the droplet can be influenced by two additional heat transfer mechanisms: heat transfer between the environment and the substrate and that between the droplet and the environment. These two are determined by the properties of the environment, including its thermal diffusivity and initial temperature. After validating the introduced LB scheme, here we examine how these parameters can affect the solidification process of a sessile droplet. We focus on two main aspects of the solidification process: the freezing time and the shape deformation of the droplet. For simulations within this section, the Stefan number is set to a constant value of 
${Ste} = 0.1$. Other parameters (shown in figure 1) are controlled through six dimensionless numbers. The parameters and corresponding ranges used in this investigation are given in table 4.
Simulation parameters and their ranges.
5.1. Freezing time
The freezing time is measured as the time 
$t$ required to fully solidify the droplet. The time is tracked through the following dimensionless time:
\begin{equation} {{Fo}} = \frac{{\alpha t}}{{l_{{{ref}}}^2}}. \end{equation} To determine which parameters from table 4 have the greatest impact on the total freezing time and how they may interact with each other, we utilize the design of experiments approach in two stages: screening and characterization. The screening step seeks to discover the most important simulation parameters while the characterization step captures the interaction effects. In the screening phase, a resolution IV standard two-level fractional-factorial design is used for six simulation parameters (listed in table 4), which consist of 16 simulations. The final shape of a solidified droplet is significantly affected by its expansion rate, but freezing time is not affected in the same way. For now, 
$\lambda$ is assumed to be zero in assessing the freezing time, and the investigation of expansion rate effects is deferred to the next section.
The statistical analysis of results are summarized in a Pareto chart in figure 6(a). The chart displays and ranks parameter effects on the freezing time. Positive effects increase the freezing time while negative effects decrease it. According to figure 6(a), the solidification time is primarily influenced by the thermal diffusivity ratio between components 
$B$ and 
$A$ (
${\alpha _{{{B(S)/A}}}}$), the contact angle (
${\psi _{{{wall}}}}$) and the initial temperature of the surrounding fluid (
${\theta _{{{ini(A)}}}}$). Other parameters fall below the 
$t$-value limit and can be considered insignificant. Although the Pareto chart includes interaction effects, the resolution for this screening phase is not high enough to clearly detect which factors are involved in these interactions.
(a) Pareto chart showing screening of simulation parameters’ effects on the freezing time. (b) Pareto chart showing characterization of four factors and their interactions to assess their effects on the freezing time.
After throwing out the unimportant factors by holding them fixed at 
${{Bo}} = 0$, 
${\alpha _{{{B(S)}}}} = {\alpha _{{{B(L)}}}} = {\alpha _{{B}}}$ and 
${\theta _{{ini(B)}}} = 0$, the study carried on with three significant parameters of 
${\alpha _{{{B/A}}}}$, 
$\theta _{{ini(A)}}$ and 
${\psi _{{{wall}}}}$ to the characterization phase. The parameter study is conducted by a two-level full factorial design with one central point of factors (
${2^3} + 1 = 9$ runs) to detect the two-factor interactions. The simulation design and corresponding freezing time are listed in the Appendix. The results displayed in figure 6(b) reveal that the only interaction that significantly affects the freezing time is between contact angle and thermal diffusivity ratio (
${\alpha _{{{B/A}}}} - {\psi _{{{wall}}}}$).
In the simulation the substrate has the lowest temperature in the computational domain. A significant portion of thermal energy of the droplet is extracted by the substrate (leading to the primary solidification). The temperature of the surrounding fluid (i.e. the environment) also changes through heat transfer with the substrate and the droplet. When the temperature of the environment drops below the melting point, a secondary solidification takes place at the interface between the droplet and the environment. The following subsections will describe how parameters selected in figure 6(b) can influence the primary and the secondary solidification and, thus, the freezing time of the droplet.
5.1.1. The thermal diffusivity ratio
When a semi-infinite liquid is placed next to a cold substrate with a constant temperature below the melting point, the thickness of the solid layer follows this equation for one-dimensional (1-D) solidification (Watanabe et al. Reference Watanabe, Kuribayashi, Honda and Kanzawa1992):
\begin{equation} {z^*} = \frac{z}{{{l_{{{ref}}}}}} = \sqrt {2 {{Ste}} {{Fo}}}. \end{equation}
The 1-D solution can still provide a good approximation for the droplet solidification when the heat exchange between droplet and environment is negligible. However, it must be noted that the environment will also be cooled down at a rate that can be different from the droplet cooling rate. This means that the surrounding fluid can reach the wall temperature much quicker than the droplet. In such situations, the freezing process is accelerated by the secondary solidification at the interface between the liquid and the environment. In this section, two different environments are used to track the solidification of a droplet. Both environments are initiated by a temperature equal to the melting point of the droplet (
${\theta _{{{ini(A)}}}} = 0$), but they are different in thermal diffusivity ratio (
${\alpha _{{{B/A}}}} = 4.0,0.25$).
The solidification of this droplet in two environments, 
${\alpha _{{{B/A}}}} = 4.0$ and 
${\alpha _{{{B/A}}}} = 0.25$, is compared in figure 7. When the thermal diffusivity of the environment is low (
${\alpha _{{{B/A}}}} = 4.0$), the surrounding fluid reacts slowly to the cold temperature. According to figure 7(a), for 
${Fo} = 1$, the temperature remained near the melting point in most of the domain, resulting in minimal secondary solidification. The solid layer has developed parallel to the substrate and is in agreement with the 1-D analytical solution. In contrast, when the rate of heat transfer inside the environment is high (
${\alpha _{{{B/A}}}} = 0.25$), the surrounding fluid cools down quickly, making secondary solidification significant. The solid interface curves upward at the edge while remaining flat at the centre of the droplet. Eventually, the concave shape of the solid layer becomes deeper until it almost forms a solid shell all over the drop at 
${Fo}=1.4$ (figure 7b). At this time more than 65 % of the droplet is solidified while in the other environment (
${\alpha _{{{B/A}}}} = 4.0$), the solid layer still has kept its flat form and the solidified percentage is about 40 %. This indicates the high dependence of freezing time on the thermal diffusivity of the environment. Higher thermal diffusivity ratio (lower thermal diffusivity in the environment) results in weaker secondary solidification and a longer freezing process (a positive effect shown in figure 6b).
Temperature distribution and solid shape of the droplet for 
${\alpha _{{{B/A}}}} = 4.0$ and 
${\alpha _{{{B/A}}}} = 0.25$ at (a) 
${Fo}=1.0$ and (b) 
${Fo}=1.4$. Other simulation parameters are set as 
${\theta _{{{ini(A)}}}} = 0$ and 
${\psi _{{{wall}}}} = 2{\rm \pi} /3$.
5.1.2. The contact angle
In this study droplets lose energy through two interfaces: one with the substrate and the other with the surrounding fluid. For ease of reference, we will refer to these two as the ‘lower’ and ‘upper’ surfaces of the droplet, respectively. Considering the droplet as a spherical cap it is easy to show that, for a given volume of the droplet, the differences in the upper surface caused by different contact angles are insignificant compared with those in the lower surface (the maximum variations of the upper and lower surface for the contact angle range of 
${\rm \pi} /3$ to 
$2{\rm \pi} /3$ are 8 % and 104 %, respectively). In the same environment, when the contact angle value is higher, the lower contact surface is reduced, resulting in a higher freezing time (a positive effect in the Pareto chart).
Figure 8 shows the variation of the freezing time as a function of contact angle for two environments with thermal diffusivity ratios of 
${\alpha _{{{B/A}}}} = 4.0$ and 
${\alpha _{{{B/A}}}} = 0.25$. In both environments the freezing time increases with higher contact angle. For this reason, hydrophobic surfaces are generally preferred as a means of delaying ice formation. But the main point is how the environment can affect the anti-ice performance of hydrophobic surfaces. When the environmental thermal diffusivity is low (
${\alpha _{{{B/A}}}} = 4.0$), changing the contact angle from 
${\rm \pi} /3$ to 
$2{\rm \pi} /3$ increases the freezing time by 155 %. However, the contact angle variation from its minimum to its maximum level only changes the freezing time by 43 % when the thermal diffusivity of the environment is high (
${\alpha _{{{B/A}}}} = 0.25$ shown as a dashed line in figure 8). This shows a strong interaction between the effects of the thermal diffusivity ratio and the contact angle (
${\alpha _{{{B/A}}}} - {\psi _{{{wall}}}}$) on the freezing time, which is highlighted before in the Pareto chart (figure 6b) as the fourth-ranked factor.
Freezing time variation under different contact angles for 
${\theta _{{{ini(A)}}}} = 0$, when the thermal diffusivity ratio is set to 
${\alpha _{{{B/A}}}} = 4.0$ (solid line) and 
${\alpha _{{{B/A}}}} = 0.25$ (dashed line).
Figure 9 compares droplet solidification on hydrophobic and hydrophilic surfaces in two environments with a similar initial temperature but different thermal diffusivity. When the environment fluid's thermal diffusivity is low (
${\alpha _{{{B/A}}}} = 4.0$ shown on the left side of figure 9), the secondary solidification is weak, and the solidified layer has a flat surface at the top for both hydrophobic and hydrophilic substrates. The difference arises when an environment fluid with a high thermal diffusivity surrounds the droplet (
${\alpha _{{{B/A}}}} = 0.25$ shown on the right side of figure 9). The secondary solidification in this case leads to a concave shape of the solidified portion of the droplet. From figure 9, it is clear that the solid interface transformation on a hydrophobic surface (upper row) is more noticeable than on a hydrophilic surface (lower row). This can be attributed to the difference in the relative strength of secondary solidification for these two droplets.
Temperature distribution and solid shape of two droplets on hydrophobic (a) and hydrophilic (b) surfaces for 
${\alpha _{{{B/A}}}} = 4.0$ (left) and 
${\alpha _{{{B/A}}}} = 0.25$ (right) when the initial temperature is set to 
${\theta _{{{ini(A)}}}} = 0$ and the dimensionless time is 
${Fo} = 1$.
When a droplet is placed on a hydrophilic surface, it tends to spread over the surface horizontally keeping a significant portion of the liquid in direct contact with the cold substrate. As a result, primary solidification dominates, even when the environment has high thermal diffusivity. In this condition, secondary solidification only has a minor impact on the liquid–solid interface near the droplet's edge. On the other hand, hydrophobic surfaces reduce primary solidification by decreasing the contact area between the liquid and the cold substrate. This exposes the droplet to a thicker layer of the environment liquid, providing more opportunities for secondary solidification to influence the freezing time. This means that using hydrophobic surfaces results in a higher sensitivity to secondary solidification. Now the question is to what extent the thermal diffusivity ratio and secondary solidification can affect the performance of hydrophobic surfaces over hydrophilic ones.
Figure 10 shows the ratio of freezing time between droplets with contact angles of 
$2{\rm \pi} /3$ and 
${\rm \pi} /3$, as a function of thermal diffusivity ratio. By increasing the thermal diffusivity of the surrounding fluid (decreasing the 
${\alpha _{{{B/A}}}}$ value), the hydrophobic surface's advantage over the hydrophilic ones decreases until it becomes almost constant above unity. It means a higher thermal diffusivity of the environment can diminish the privilege of a hydrophobic surface, but it cannot completely remove it. Because in this set-up where the substrate has the lowest temperature, the secondary solidification is always weaker than what is caused by the cold wall. When the solidification from the lower surface is dominant, a higher contact angle will always result in a longer freezing time.
Variation of the freezing time ratio between droplets with contact angles of 
$2{\rm \pi} /3$ and 
${\rm \pi} /3$ as a function of the thermal diffusivity ratio when the initial temperature is 
${\theta _{{{ini(A)}}}} = 0$.
5.1.3. The initial temperature of the environment
When the environment temperature is lower than the melting point, secondary solidification speeds up the freezing process. Conversely, in an environment with high temperatures, freezing time will be longer (a positive effect shown in the Pareto chart, figure 6b). The relative strength of secondary solidification determines the importance of the initial temperature of the environment in freezing time. In the current set-up, solidification primarily occurs at the bottom of the droplet. Therefore, the initial temperature of the environment is not expected to be ranked in the top tier of the Pareto chart. Conditions can change in a different set-up.
When the temperature difference between the droplet and its surroundings is significant enough to disregard solidification from the substrate, the primary factor affecting the freezing time is the environment. The upper surface of the droplet is where this influence is felt the most. In order to simulate such a condition, the substrate is assumed to be thermally isolated. The environment is initiated at 
${\theta _{{{ini(A)}}}} = - 1$, the temperature at the upper boundary of the domain is kept constant at 
$\theta = - 1$, and three droplets with wall contact angles of 
${\rm \pi} /3$, 
${\rm \pi} /2$ and 
$2{\rm \pi} /3$ are rested at the substrate to freeze by the environment. The solid shape of the three droplets at times 
${Fo} = 2$ and 
${Fo} = 4$ can be seen in figure 11. Because the substrate is treated as an adiabatic boundary, the inner surface of the solid layer of the droplet is perpendicular to the lower surface (substrate). The figure shows that the solidification process occurs at equal rates in all droplets. That is because the variation of the upper surface due to different contact angles is not significant. Finally, three droplets solidify at 
${{{Fo}}_{{{freezing(}}{\rm \pi} {{/3)}}}} = 4.81$, 
${{{Fo}}_{{{freezing(}}{\rm \pi} {{/2)}}}} = 4.48$ and 
${{{Fo}}_{{{freezing(2}}{\rm \pi} {{/3)}}}} = 4.28$. It means that when solidification occurs primarily on the upper surface, higher contact angles lose effectiveness in anti-icing and can even decrease the freezing time. In this simulation, the condition is achieved by isolating the substrate. However, other phenomena, such as surface evaporation or wind, can make solidification at the upper surface dominant and lead to the same results.
The evolution of solid layers for three droplets with contact angles of 
${\rm \pi} /3$, 
${\rm \pi} /2$ and 
$2{\rm \pi} /3$ when the solidification initiated from the upper surface. Results are shown for (a) 
$Fo=2.0$. (b) 
$Fo=4.0$.
5.2. Shape deformation
During solidification, the volume of a droplet changes if its expansion rate is not zero. The final shape of the solidified droplet is determined by the evolution of the solid–liquid interface. At the early stage of solidification, a flat layer of solid grows next to the wall. In this stage, the heat transfer is limited to the cold substrate and the interface stays fairly flat. As the solid layer grows, the shape of the liquid–solid interface becomes sensitive to edge conditions where the heat transfer is conducted through a tri-junction between the solid, liquid and surrounding fluid (figure 12). Here we track the solid front through variation of two heights: the solid thickness at the centre of the droplet (
${h_c}$) and the height of the solid front at the tri-junction (
${h_e}$), both are shown in figure 12. The two heights can be made dimensionless by dividing them by the initial height of the droplet, denoted as 
$H$.
Schematic of a sessile droplet during its solidification.
Section 5.1 showed how the properties of the environment affect the secondary solidification at the tri-junction. Thermal diffusivity was found to be the most significant factor in the solidification process. Here, our discussion focuses on how the thermal diffusivity ratio (in the range of 
$0.5$–
$20$) can affect the final shape of a solidified droplet. For this part, the initial temperature of the environment is set to 
${\theta _{{{ini(A)}}}} = 0$, and the expansion rate is fixed at 
$\lambda = 10\,\%$.
Figure 13(a) shows the evolution of 
${h_c}$, and 
${h_e}$ as functions of the freezing time fraction (
$\tau = {{Fo}}/{{F}}{{{o}}_{{{freezing}}}}$) for 
${\alpha _{{{B/A}}}} = 20$. Also included in this figure are the component index distribution and droplet shapes before and after solidification. This figure indicates that the solid layer remains flat until 
$\tau = 0.3$. After that, 
${h_e}$ tends to be slightly higher than 
${h_c}$, but stays close to 
${h_c}$ for much of the solidifying process. Finally, the solidification is finished with 10 % extra volume. There is no sharp tip above the droplet, but concentration of expanded volumes on the central part of the droplet leads to a conical shape at the top. The final shape of the droplet is influenced by the gradual expansion of volume. For a period of time, the solid layer grows in a flat shape, allowing the extra volume induced by solidification to expand horizontally. This ensures smooth shape variation without severe shape deformation at the top of the droplet. The ice–liquid interface deformation from a flat to a concave shape can occur sooner by increasing the thermal diffusivity of the environment. At 
${\alpha _{{{B/A}}}} = 11.5$, the solid layer stays flat for just 21 % of the solidification time. The concave shape of the solid interface holds the extra volume in the central part and ultimately leads to a sharper tip above the droplet.
Evolution of the solid layer as a function of freezing time fraction, and the component index distribution in the solidified droplet for (a) 
${\alpha _{{{B/A}}}} = 20$ and (b) 
${\alpha _{{{B/A}}}} = 11.5$ when the expansion rate is 
$\lambda = + 10\,\%$.
5.2.1. Solid shell formation and residual thermal stress
As the thermal diffusivity of the environment increases, the edge of the droplet experiences faster solid growth. Eventually, a solid shell develops and prevents any further deformation during solidification. Figure 14(a) shows such a condition where the thermal diffusivity ratio is set to 
${\alpha _{{{B/A}}}} = 1$. According to this figure, the solid shell is formed at 
$\tau =0.8$, while the central height is 
${h_c} = 0.69H$. In this stage, shape deformation stops, and the core liquid's tendency to change volume may lead to two outcomes. The solid shell may be strong enough to prohibit the core liquid from further volume change, leading to residual stress inside the droplet. This is the case in our simulation in figure 14(a). At 
$\tau =1$ when the solidification completes, the drop has a sharp tip but its volume is only 8 % larger than its initial shape. Forming a solid shell has another possibility for the final shape. Solidification of liquid inside the shell may induce cracks in the outer solid layers. Some liquid can pass through the cracks to solidify on the surface. The final shape will be slightly warped, probably with a tip at the top. This kind of deformation cannot be modelled using the present scheme, because the solid shell in simulations is rigid and unbreakable.
Evolution of the solid layer as a function of freezing time fraction, and the component index distribution in the solidified droplet for (a) 
${\alpha _{{{B/A}}}} = 1$ and (b) 
${\alpha _{{{B/A}}}} = 0.5$ when the expansion rate is 
$\lambda = + 10\,\%$.
When the solid outer shell is strong enough to prevent any further volume change, there is a chance of residual stresses after solidification. This is a very common situation that can happen in manufacturing processes. The cooling rate of a casting product is not uniform. Usually, the temperature of surface layers decreases more rapidly than the core temperature. When the expansion coefficient of the molten material is not zero, any expansion or contraction of the core part during the solidification will be prevented by the solid outer layers resulting in localized residual compressive or tensile stress. The calculation of the value and location of these stresses are rather complex depending on the geometry, temperature condition and mechanical properties of the component. Full analysis of these stresses can be achieved by experimental data and elastoplastic analysis, which is outside the scope of this study. However, the introduced LB model has a specific feature that can be used to predict these residual thermal stresses.
Section 2 introduces a new source term for the convective Cahn–Hilliard equation that simulates volume change by manipulating the component index around its equilibrium value. The increased or decreased component index will then be balanced by moving components into or out of the solidification interface. If the deviated component index value is surrounded by a solidified domain, the solid treatment introduced in § 3.3 prevents any further component movement. This means that the deviated component index value remains trapped inside the solid domain during the simulation. Emphasizing that residual thermal stress is caused by blockage of volume change of a core part by an outer solid layer, we can use residually unbalanced component index distribution as a sign of possible residual thermal stresses. Note that in the case of severe contraction, the method must be used with caution. Because the unbalanced component index inside a solid shell may fall below the zero value that is not physically valid.
In figures 13 and 14 the final shape of the droplets is also represented by the distribution of the component index. For a thermal diffusivity ratio of 
${\alpha _{{{B/A}}}} = 1$, the solid shell forms at 
$\tau =0.8$. The shape deformation stops and the solidification of the remaining liquid inside the droplet results in high index values (depicted as red regions) that resembles the residual compressive stresses. Note that there are no such red regions in figure 13.
Despite full solidification, the droplet can still maintain its spherical cap shape if the thermal diffusivity ratio is sufficiently low. Figure 14(b) illustrates the process of solidifying a droplet when the thermal diffusivity ratio is 
${\alpha _{{{B/A}}}} = 0.5$. The solid shell forms at 
$\tau = 0.65$ when the volume is expanded by only 6 %. The drop completes its solidification with this percentage and a shape that has retained its cap shape. This condition has a larger red zone than figure 14(a), which suggests a greater chance of residual thermal stresses.
5.2.2. Sign of volume change
So far, we have described three possible shape deformations for a droplet that expands during solidification. A ‘uniform’ expansion occurs when the environment has a low thermal diffusivity (figure 13a). A ‘sharp’ deformation happens when the solid–liquid interface forms a concave shape and maintains this shape for most of the freezing process (figures 13b and 14a). An ‘incomplete’ deformation is caused by the creation of a solid shell in the presence of a high environmental thermal diffusivity (figure 14b). A droplet that contracts during solidification can experience any of these possible deformations. However, with the present scheme, only uniform and sharp deformations can be modelled. As previously mentioned, an incomplete deformation may result in physically invalid results when there is severe contraction.
Figure 15 illustrates the solid layer's evolution during a uniform deformation while a droplet with a 
$\lambda =-10\,\%$ is surrounded by an environment with a thermal diffusivity ratio of 
${\alpha _{{{B/A}}}} = 8$. A flat solid layer grows until 
$\tau = 0.18$. The interface maintains its flat shape even after the solid edge starts to grow faster (
$\tau =0.4$). At the end, the volume of solidified droplet is reduced by 
$10\,\%$, but in a uniform manner. A sharp deformation can be achieved with a higher thermal diffusivity of the environment. Such a condition is shown in figure 16 where the thermal diffusivity ratio is 
${\alpha _{{{B/A}}}} = 0.5$. As shown in this figure, the liquid at the edges solidifies faster than the central region during the early stages of solidification. At 
$\tau = 0.4$, the interface takes on a concave shape, as opposed to the semi-flat interface seen in figure 15. Finally, when the solidification is completed, a dimple forms at the top of the droplet.
Evolution of the solid layer as a function of freezing time fraction, and the component index distribution in the solidified droplet for 
${\alpha _{{{B/A}}}} = 8$ when the expansion rate is 
$\lambda = - 10\,\%$.
Evolution of the solid layer as a function of freezing time fraction, and the component index distribution in the solidified droplet for 
${\alpha _{{{B/A}}}} = 0.5$ when the expansion rate is 
$\lambda = - 10\,\%$.
6. Conclusion
In this paper we introduced a new LB model to simulate solidification processes in a binary multicomponent system. Based on the free energy approach, the model uses a modified version of the convective Cahn–Hilliard equation to capture the density variation as well as changes in other properties during solidification. The model was validated against analytical, numerical and experimental results available in the literature. The results demonstrated that the new LB framework is capable of simulating liquid–solid phase transitions with volume expansion or contraction due to changes in density.
Using the new model, the solidification of a sessile droplet was investigated in detail. Six different parameters are examined to see how they – individually or combined – can affect the freezing time of the droplet. The results showed that the environmental conditions and solid surface properties have the most dominant influence on the freezing time. By analysing the interaction effects of these parameters, we found that the effects of surface hydrophobicity can be partially diminished or intensified by the thermal diffusivity of the environment. It is also shown that the freezing times are not significantly different on hydrophobic and hydrophilic surfaces when solidification primarily occurs at the droplet–environment interface.
The variation of droplet shape was also simulated using the new model. In the case of volume expansion, numerical results indicated that the droplet tends to develop a cone on the upper surface. Conversely, volume contraction can lead to a dimple at the top of the solidified droplet. The new LB model also provides a new mechanism to assess the possibility of thermal residual stresses that may occur during solidification. The importance of environmental conditions was emphasized once again in the simulation of droplet solidification, taking into account both volume expansion and contraction. The results obtained from this work may be of great significance to better understand the solidification mechanisms when the liquid–solid phase transition occurs among a surrounding fluid.
Acknowledgements
Thanks to Dr K. Shi for providing experimental data. The authors would like to acknowledge his contribution.
Declaration of interests
The authors report no conflict of interest.
Appendix. Simulation designs and results obtained from simulations
This appendix lists simulation designs and results obtained from simulations.
Simulation design and freezing time obtained from simulations in the screening phase.
Freezing time obtained from simulations for the characterization phase.
