A consistent and conservative Phase-Field model for thermo-gas-liquid-solid flows including liquid-solid phase change

02/13/2021
by   Ziyang Huang, et al.
Purdue University
0

In the present study, a consistent and conservative Phase-Field model is developed to study thermo-gas-liquid-solid flows with liquid-solid phase change. The proposed model is derived with the help of the consistency conditions and exactly reduces to the consistent and conservative Phase-Field method for incompressible two-phase flows, the fictitious domain Brinkman penalization (FD/BP) method for fluid-structure interactions, and the Phase-Field model of solidification of pure material. It honors the mass conservation, defines the volume fractions of individual phases unambiguously, and therefore captures the volume change due to phase change. The momentum is conserved when the solid phase is absent, but it changes when the solid phase appears due to the no-slip condition at the solid boundary. The proposed model also conserves the energy, preserves the temperature equilibrium, and is Galilean invariant. A novel continuous surface tension force to confine its contribution at the gas-liquid interface and a drag force modified from the Carman-Kozeny equation to reduce solid velocity to zero are proposed. The issue of initiating phase change in the original Phase-Field model of solidification is addressed by physically modifying the interpolation function. The corresponding consistent scheme is developed to solve the model, and the numerical results agree well with the analytical solutions and the existing experimental and numerical data. Two challenging problems having a wide range of material properties and complex dynamics are conducted to demonstrate the capability of the proposed model.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 8

page 17

page 18

page 19

page 22

page 23

01/12/2021

A consistent and conservative model and its scheme for N-phase-M-component incompressible flows

In the present work, we propose a consistent and conservative model for ...
03/14/2021

Implementing contact angle boundary conditions for second-order Phase-Field models of wall-bounded multiphase flows

In the present work, a general formulation is proposed to implement the ...
10/15/2021

A fully conservative sharp-interface method for compressible mulitphase flows with phase change

A fully conservative sharp-interface method is developed for multiphase ...
07/27/2020

A staggered-projection Godunov-type method for the Baer-Nunziato two-phase model

When describing the deflagration-to-detonation transition in solid granu...
12/26/2019

Travelling wave mathematical analysis and efficient numerical resolution for a one-dimensional model of solid propellant combustion

We investigate a model of solid propellant combustion involving surface ...
03/19/2022

Meshfree One-Fluid Modelling of Liquid-Vapor Phase Transitions

We introduce a meshfree collocation framework to model the phase change ...
05/05/2020

Direct numerical simulation of the combustion of a suspended droplet in normal gravity

DropletSMOKE++ is a multiphase CFD framework based on OpenFOAM, original...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Liquid-Solid phase change (or solidification/melting) and its interaction with the surrounding air are ubiquitous in various natural and/or industrial processes, e.g., latent thermal energy storage (LTES) systems (Shmuelietal2010; VogelThess2019), welding (Pitschenederetal1996; Chanetal1984; Zhaoetal2011; Saldi2012), casting (DantzigRappaz2016; Huangetal2018), and additive manufacturing (AM) (Panwisawasetal2017; Heetal2020; Linetal2020). This motivates researchers to develop physical and high-fidelity models to further understand the complex physics and dynamics, and accurately predict the behavior of materials in order to produce high-quality products. Such a problem includes many challenging factors, such as a wide range of material properties, evolution of the liquid-solid interface due to solidification/melting, deformation of the gas-liquid interface due to fluid (gas/liquid) motions and surface tension, heat transfer that drives the phase change, and fluid-structure interaction between the fluid and solid, and all these factors are coupled and can be influential. In the present study, we call the problem thermo-gas-liquid-solid flows with liquid-solid phase change.

In spite of the complexity of the problem, it can be separated into two basic problems which are the two-phase incompressible flow and the solidification with convection. Several studies focused on these individual problems. For two-phase incompressible flow, the front-tracking method UnverdiTryggvason1992; Tryggvasonetal2001, the level-set method OsherSethian1988; Sussmanetal1994; SethianSmereka2003; Gibouetal2018, the conservative level-set method OlssonKreiss2005; Olssonetal2007; ChiodiDesjardins2017, the volume-of-fluid (VOF) method HirtNichols1981; ScardovelliZaleski1999; OwkesDesjardins2017, the THINC method Xiaoetal2005; Iietal2012; XieXiao2017; Qianetal2018, and the Phase-Field (or Diffuse-Interface) method Andersonetal1998; Jacqmin1999; Shen2011; Huangetal2020 have been developed to locate different phases. The smoothed surface stress method Gueyffieretal1999, the continuous surface force (CSF) Brackbilletal1992, the ghost fluid method (GFM) Fedkiwetal1999; Lalanneetal2015, the conservative and well-balanced surface tension model Abu-Al-Saud2018, and the Phase-Field method derived from the energy balance or the least-action principle Jacqmin1999; Yueetal2004 have been developed to model the surface tension, and a balanced-force method (Francoisetal2006) is proposed to incorporate the surface tension model to the fluid motion. Interested readers should refer to ProsperettiTryggvason2007; Tryggvasonetal2011; Shen2011; Mirjalilietal2017; Popinet2018. For the solidification problem with convection, one of the most popular methods is the enthalpy-porosity technique (Volleretal1987; VollerPrakash1987; Brentetal1988; VollerSwaninathan1991; RoslerBruggemann2011), which can be implemented in a fixed grid. The liquid fraction of the phase change material is algebraically related to the local temperature. Therefore, the enthalpy change due to phase change can be evaluated and becomes a source in the energy equation. A drag force proportional to the velocity is added to the momentum equation to stop the solid motion, which is the same idea as the fictitious domain Brinkman penalization (FD/BP) method for the fluid-structure interaction (Angotetal1999; BergmannIollo2011). The most popularly used drag force in the enthalpy-porosity technique is modified from the Carman-Kozeny equation (Carman1997) for the porous medium. Another popular method for solidification is the Phase-Field method (Boettingeretal2002; Chen2002; Echebarriaetal2004; KimKim2005; Tanetal2011; Jietal2018; Luetal2018), where the liquid-solid interface has a small but finite thickness. Different from the enthalpy-porosity technique using an algebraic relation, the Phase-Field method introduces an additional equation to govern the evolution of the liquid fraction of the phase change material, and therefore is flexible to include more complicated physics, e.g., anisotropy and solute transport in an alloy. After coupling the Phase-Field method with the hydrodynamics, the melt convection can be modeled (Nestleretal2000; Beckermannetal1999; ChenYang2019; ZhangYang2020). Other methods for modeling the liquid-solid phase change are reviewed in (SalcudeanAbdullah1988; Samarskiietal1993; Voller1996; HuArgyropoulos1996; Dutiletal2011; DhaidanKhodadadi2015; Sultanaetal2018).

Most existing models for the thermo-gas-liquid-solid flows with liquid-solid phase change follow a similar procedure: The deforming gas-liquid interface is located by a certain interface capturing method and the continuous surface tension force is added, while the enthalpy-porosity technique is directly applied without any further changes to adapt to the appearance of a new phase. The volume-of-fluid method is the most popular choice and is used, e.g., in (Shmuelietal2010; Saldi2012; Kimetal2013; Yanetal2017; Panwisawasetal2017; VogelThess2019; Heetal2020). The level-set and conservative level set methods are recently used in (Yanetal2018) and (Linetal2020), respectively. Some additional physics are introduced to the models, e.g., the thermo-capillary effect (Panwisawasetal2017; Yanetal2018; Heetal2020; Linetal2020) and recoil pressure (Panwisawasetal2017; Heetal2020; Linetal2020). Another recent model (Zhangetal2020) follows the same strategy but instead uses the Phase-Field model in (Ramirezetal2004) for anisotropic solidification and the conservative Phase-Field method (ChiuLin2011) as the interface capturing method. In spite of its widespread applications, such a well-accepted modeling strategy has the following critical issues. (i) The volume fractions of the phases are ambiguously defined. In the models applying the enthalpy-porosity technique, the liquid fraction is meaningful only inside the phase change material, while it directly appears in the energy equation defined in the entire domain including the gas phase. As a result, the meaningless value of the liquid fraction in the gas phase is also counted in the energy equation. Another example is in (Zhangetal2020), where two liquid fractions are defined for the same liquid phase, one from the solidification model and the other from the interface capturing method. Since the solidification model and the interface capturing method have no explicit/direct connection, these two liquid fractions may inconsistently label the liquid location. (ii) The surface tension and drag forces can appear at the wrong location because the volume fractions, which are not clearly defined, are needed to compute the forces. Based on the formulations, e.g., in (Saldi2012; Panwisawasetal2017; Yanetal2018; Linetal2020; Heetal2020), the surface tension at the gas-liquid interface will falsely appear at the gas-solid interface, and the drag force will falsely appear in the gas phase when the local temperature is lower than the solidus temperature. Some artificial operations need to be added but they have seldom been mentioned in the literature. An exception is in (Zhangetal2020) where a bounce-back scheme near the gas-solid interface is employed since the gas-solid interface is unable to be effectively labeled by the model, but details of the bounce-back scheme are not provided. (iii) Physical principles can be violated, depending on the material properties. The most obvious example is the mass conservation. In, e.g., (RoslerBruggemann2011; Panwisawasetal2017; Yanetal2017; Yanetal2018; Zhangetal2020; Heetal2020; Kimetal2013), the velocity is divergence-free, implying that both the volumes of the gas and phase change material will not change. This restricts the application of the model to problems having matched liquid and solid densities. However, the problems studied in (Kimetal2013; Linetal2020) are outside that category, and therefore the volume of the phase change material needs to change in order to satisfy the mass conservation. It should be noted that the volume change in (Linetal2020) is from evaporation, not solidification/melting, and the velocity is divergence-free without evaporation. As a result, the divergence-free velocity is contradicting the mass conservation. Although the studies in (Shatikianetal2005; Shmuelietal2010; Hosseinizadehetal2011; VogelThess2019) captures the volume change, detail formulations, i.e., the divergence of the velocity, are not provided. Other physical principles, e.g., the momentum and energy conservation and the Galilean invariance of the models, have never been examined. (iv) The problem of interest often has a large density ratio, which can be between the gas and liquid. It has been well-known that the so-called consistent method Rudman1998; Bussmannetal2002; ChenadecPitsch2013; OwkesDesjardins2017; RaessiPitsch2012; Nangiaetal2019; Xieetal2020; Huangetal2020; Huangetal2020CAC needs to be implemented to produce physical results for multiphase flows, while this has never been considered in the existing models for the problem of interest.

The aforementioned issues in the existing models for the problem of interest are originated in simply “combining”, not physically “coupling”, the models for solidification and two-phase flow. In the present study, those critical issues are properly addressed and a consistent and conservative Phase-Field model is developed for the thermo-gas-liquid-solid flows with liquid-solid phase change. All the dependent variables are defined in a fixed regular domain, which is convenient for numerical implementation.

The novelty of the present study is multi-fold and the proposed model enjoys the following physical properties:

  • In deriving the proposed model, several consistency conditions proposed in our previous works (Huangetal2020; Huangetal2020CAC; Huangetal2020N; Huangetal2020B; Huangetal2020NPMC) are considered. Our earlier works do not include phase change or temperature variations in the fluids, and the present study is the first implementation of the consistency conditions to phase change problems, which further demonstrates their generality.

  • The proposed model exactly recovers the consistent and conservative Phase-Field method for incompressible two-phase flows (Huangetal2020) when the solid phase is absent, the fictitious domain Brinkman penalization (FD/BP) method for fluid-structure interactions (Angotetal1999; BergmannIollo2011) when the liquid-phase is absent, and the Phase-Field model of solidification of a pure material (Boettingeretal2002) when the gas phase is absent.

  • The proposed model defines the volume fractions of individual phases unambiguously and ensures their summation to be unity everywhere. The local mass conservation is strictly satisfied, from which the divergence of the velocity is non-zero, and therefore the volume change during solidification/melting is captured.

  • The momentum is conserved when the solid phase is absent, and the no-slip condition at the solid boundary results in the momentum change. The energy conservation and Galilean invariance are also satisfied by the proposed model.

  • The momentum transport is consistent with the mass transport of the gas-liquid-solid mixture, which greatly improves the robustness of the model for large-density-ratio problems and avoids unrealistic interface deformation.

  • Isothermal (or temperature equilibrium) solutions are admissible, thanks to satisfying the consistency conditions, which prevents producing any fictitious fluctuations of the temperature.

  • Novel continuous surface tension and drag force models are proposed, which are activated only at proper locations. The interpolation function in the solidification model (Boettingeretal2002) is modified to include the capability of initiating phase change when there is only the liquid/solid-state of the phase change material.

These physical properties of the proposed model are independent of the material properties.The proposed model is validated and its capability is demonstrated using a consistent numerical scheme. The proposed model includes all the basic ingredients and challenging aspects of the problem, and additional physics can be incorporated conveniently following the same framework.

The rest of the paper is organized as follows. In Section 2, the proposed model and its properties are elaborated in detail. In Section 3, the numerical procedure to solve the proposed model is introduced. In Section 4, various numerical tests are performed to validate the properties of the proposed model, and challenging problems are simulated to demonstrate the capability of the proposed model. In Section 5, the present study is concluded.

2 Governing equations

The problem considered includes two materials, which are a gas “” and a phase change material “” experiencing solidification or melting. Therefore, in the entire domain , there are three phases: the gas phase including only “”, and the liquid and solid phases of “”. The liquid-solid phase change is driven by temperature. The part of occupied by “” is denoted by , and similarly, we use , , and to denote the domains occupied by “”, the liquid phase, and the solid phase of “”, respectively. As a result, we have . In addition, boundaries of the domains are denoted with “” in front of the corresponding domains, for example, is the boundary of . The material properties of the gas phase “” and the liquid and solid phases of “” are assumed to be constant and denoted by , , and , respectively, where can be density , viscosity , specific heat , and heat conductivity .

The proposed consistent and conservative model for thermo-gas-liquid-solid flows including liquid-solid phase change is elaborated in Section 2.1. Then in Section 2.2, the physical properties and the relations of the proposed model to some other multiphase models are analyzed. When the proposed model is derived in Section 2.1, several consistency conditions will be applied so that the proposed model is able to produce physical results. The definitions of the consistency conditions are:

  • Consistency of reduction: The multiphase system should be able to recover the corresponding systems including fewer phases.

  • Consistency of volume fraction conservation: The phase change equation should be consistent with the governing equation for the volume fraction of “”, when “” is in a fully liquid/solid-state.

  • Consistency of mass conservation: The mass conservation equation should be consistent with the governing equation for the volume fraction of “”, the phase change equation, and the density of the multiphase mixture. The mass flux and the divergence of the velocity in the mass conservation equation should lead to a zero mass source.

  • Consistency of mass and momentum transport:

    The momentum flux in the momentum equation should be computed as a tensor product between the mass flux and the velocity, where the mass flux should be identical to the one in the mass conservation equation.

These consistency conditions are generalized from their correspondences for isothermal multimaterial incompressible flows (Huangetal2020; Huangetal2020CAC; Huangetal2020N; Huangetal2020B; Huangetal2020NPMC) to include effects of phase change and temperature variation.

It should be noted that any thermo-gas-liquid-solid problems locally are gas-liquid, gas-solid, or liquid-solid problems and can be isothermal. Therefore, all these circumstances should be admissible by a thermo-gas-liquid-solid model in order to produce correct physical dynamics, which is emphasized by the consistency of reduction. To physically understand the meaning of the rest of the consistency conditions, one needs to first realize that the Phase-Field method includes relative motions of the materials, modeled as a diffusive process. The consistency conditions are general principles to incorporate the mass, momentum, and energy transports due to the relative motions of the materials into the conservation equations.

2.1 The proposed model

2.1.1 The Cahn-Hilliard equation

The interfacial dynamics of “” and “” is modeled by the Cahn-Hilliard equation (CahnHilliard1958) with convection:

(1)

Here, is the order parameter of the Cahn-Hilliard equation and is considered as the volume fraction of “” in . is the velocity, whose divergence can be nonzero. and are the mobility and chemical potential of , respectively. is the mixing energy density of , which is proportional to both the thickness of “-” interface and the surface tension at the gas-liquid interface . is the double-well potential, having minimums at or , and is the derivative of with respect to . is the Phase-Field flux of , including both the convection and diffusion fluxes in the Cahn-Hilliard equation. Unless otherwise specified, the homogeneous Neumann boundary condition is applied.

The Cahn-Hilliard equation Eq.(1) is derived from the Ginzburg-Landau free energy functional using the gradient flow, and the chemical potential in Eq.(1) is the functional derivative of the Ginzburg-Landau free energy functional with respect to the order parameter . The Cahn-Hilliard equation has been widely used in modeling two-phase incompressible flows, e.g., in (Jacqmin1999; Dingetal2007; Abelsetal2012; Huangetal2020), and its derivation has already been given in many references, e.g., in (Shen2011; Fengetal2005; Yueetal2004).

2.1.2 The phase change equation

In the present study, we consider the Allen-Cahn Phase-Field model for the solidification of a pure material (Boettingeretal2002) and the convection term is added:

(2)

Notice that Eq.(2) is defined in . Here, is the order parameter of the Allen-Cahn dynamics, and is considered as the volume fraction of the liquid phase of “” in . is the mobility of . and are the Gibbs-Thomson and linear kinetic coefficients, respectively, in the Gibbs-Thomson equation of the liquid-solid interface. is the mixing energy density of , which is related to the thickness of the liquid-solid interface . is the latent heat of the liquid-solid phase change. is the temperature and is the melting temperature of the phase change. is the interpolation function, monotonically increasing from to and having extreme points at and , and is the derivative of with respect to . is the double-well potential function defined identically to the one in Eq.(1). Unless otherwise specified, the homogeneous Neumann boundary condition is applied.

This solidification/melting model Eq.(2) is the basis of many more complicated models including, e.g., components and/or anisotropy (Boettingeretal2002; Chen2002; Jietal2018; KimKim2005; ChenYang2019; ZhangYang2020). The model can be derived either thermodynamically from a free energy functional using the gradient flow or geometrically from the Gibbs-Thomson equation, and details are available in (Boettingeretal2002; Beckermannetal1999; AllenCahn1979). The first two terms in the parentheses on the right-hand side of Eq.(2) models the curvature driven effect on the phase change, while they are, at the same time, competing with each other to maintain the thickness of the liquid-solid interface. The effect of the temperature is modeled by the last term in the parentheses.

Solving Eq.(2) is very challenging because it is defined in which is evolving with time. It would be more convenient if we can obtain an equivalent equation to Eq.(2) but defined in . Therefore, the diffuse domain approach (Lietal2009) is applied, and we use , the volume fraction of “”, as an approximation to the indicator function of whose value is in but elsewhere. The equivalence of Eq.(2) in is

(3)

We directly apply the formulations of the diffuse domain approach in (Lietal2009) to obtain Eq.(3) from Eq.(2), and details of the approach are available in (Lietal2009). Two modifications will be performed to Eq.(3) to address the following two issues.

The first issue is about initiating the phase change. It should be noted that the terms in the parentheses on the right-hand side of Eq.(2), which are the driving forces for the phase change, are nonzero only at . In other words, given “” in fully solid (liquid) state at the beginning, melting (solidification) will never happen no matter how high (low) the temperature is. This issue is originated in the definition of in Eq.(2) whose extreme points are and , i.e., . These extreme points are the same as the equilibrium states of . Defining in this way, as mentioned in (Boettingeretal2002), is an improvement from using , e.g., in (Javierreetal2006), in the sense that the equilibrium state of is always or , independent of the temperature. However, defining preserves the driving force from the temperature when or . The equilibrium states of should depend on the temperature. If the temperature is larger than the melting point, the equilibrium state should be (liquid), while it should be (solid) if . To achieve this property, we propose defined in Eq.(4) in the present study, which combines the advantages of in (Boettingeretal2002) and in (Javierreetal2006), but avoids their disadvantages. Only when (liquid) and (undercool) or when (solid) and (overheat), is , so that the effect of the temperature on the phase change is included. In other cases, is the same as with . As a result, the equilibrium state of when is (solid), and it is (liquid) when .

The second issue is about the existence of fully liquid/solid state of “”. For example, given, and , we obtain from Eq.(2), which implies . In other words, Eq.(2) admits the existence of fully liquid state of “”, when the temperature is larger than the melting temperature and there is no solid phase at the beginning. This property should be inherited by Eq.(3), while this is not the case. Using the same condition, we are unable to obtain from Eq.(3). To address this issue, we apply the consistency of volume fraction conservation. After comparing Eq.(3) along with and to Eq.(1), we discover that the convection velocity in Eq.(3) needs to be replaced with . As a result, we have on the left-hand side and on the right-hand side, given and , and these two sides are equal to each other from the Cahn-Hilliard equation Eq.(1). Therefore, the existence of fully liquid state of “” is admitted by Eq.(3) after the modification. The same is also true for the solid state of “”.

After applying the above mentioned modifications to Eq.(3) to address these two issues, we obtain the phase change equation:

(4)

Here, is called the chemical potential of , is modified from to initiate the phase change, and is the Phase-Field flux of .

2.1.3 The volume fractions and material properties

Based on the Cahn-Hilliard equation Eq.(1) and the phase change equation Eq.(4), the volume fractions of the gas, liquid, and solid phases in the proposed model are defined unambiguously. The volume fraction of the gas phase (or “”) in is , and they are and for the liquid and solid phases, respectively. The volume fraction of “” in is . It is clear that is always true in the entire domain.

With the volume fractions of the phases in hand, the material properties of the gas-liquid-solid mixture and their fluxes are computed as

(5)

Here, represents a certain material property, e.g., density , and and are the Phase-Field fluxes defined in Eq.(1) and Eq.(4), respectively.

2.1.4 The mass conservation

To determine the mass transport of the model, the consistency of mass conservation is applied. First, the density of the multiphase mixture is obtained following Eq.(5)

(6)

After combining Eq.(6) with the Cahn-Hilliard equation, Eq.(1), and the phase change equation, Eq.(4), the mass of the multiphase mixture is governed by

(7)

where is the consistent mass flux defined in Eq.(5). Finally, to obtain a zero mass source in Eq.(7), the divergence of the velocity should satisfy

(8)

Eq.(8) illustrates the volume change due to the liquid-solid phase change. As a result, the mass conservation equation of the model is

(9)

Therefore, the mass of the gas-liquid-solid mixture is locally conserved, even though the phase change happens. This is achieved at the expense of changing the volume of “”, as indicated in Eq.(8). On the other hand, if the densities of the liquid and solid phases of “” are the same, the volume of “” remains the same, and therefore the divergence of the velocity is zero, which can also be seen in Eq.(8). It should be noted that the mass conservation equation Eq.(9) is not an independent equation in the model. Instead, it is derived from Eq.(5) and Eq.(8), as elaborated in this section. Therefore, we don’t need to explicitly solve the mass conservation equation Eq.(9).

2.1.5 The momentum equation

The motion of the phases (or materials) is governed by the momentum equation:

(10)

Here, is the pressure, and is the gravity. As all the gas-liquid, gas-solid, and liquid-solid interfaces are immersed in , their effects are modeled as volumetric forces, i.e., and , in the momentum equation Eq.(10). is the surface tension force, modeling the surface tension at the gas-liquid interface, while is the drag force, modeling the no-slip boundary condition at the gas-solid and liquid-solid interfaces by enforcing in the solid phase. is the drag coefficient of , and and are model parameters of . is the given solid velocity and is set to be zero unless otherwise specified. Note that the momentum is transported with the consistent mass flux that also appears in the mass conservation equation Eq.(9). This follows the consistency of mass and momentum transport, and is essential to obtain the kinetic energy conservation (when only the pressure gradient is present on the right-hand side of Eq.(10)) and the Galilean invariance, as analyzed in (Huangetal2020).

To model the surface tension at the gas-liquid interface, we apply the commonly used Phase-Field formulation , which can be derived from either the energy balance Jacqmin1999; Huangetal2020CAC or the least-action principle Yueetal2004; Shen2011. However, is activated at all “-” interfaces, including both the gas-liquid and gas-solid interfaces. To remove its contribution at the gas-solid interface, is multiplied by so that the surface tension force remains to be at the gas-liquid interface and smoothly reduces to zero away from it. As a result, we obtain in Eq.(10). To more clearly illustrate the distribution of the surface tension force, we consider a gas bubble at the center of a unit domain with a radius , the bottom half of which is in contact with the solid phase while the upper half is in contact with the liquid phase. Fig.1 schematically shows the magnitude of the surface tension forces. We will further investigate the surface tension force in Section 4.1.5. It should be noted that the thermo-capillary effect has not been considered in the present study and the surface tension is treated as a constant.


Figure 1: Magnitude of the surface tension force . a) The proposed surface tension force , b) The original surface tension force . The proposed surface tension force shown in a) removes its contribution at the gas-solid interface, compared to the original formulation shown in b).

To model the no-slip boundary condition at both the gas-solid and liquid-solid interfaces, the velocity inside the solid phase needs to be the given value . We add a drag force, formulated as , to the momentum equation Eq.(10), where the drag coefficient can be considered as the Lagrange multiplier that enforces inside the solid phase. should depend on the volume fraction of the solid phase in such a way that it is zero away from the solid phase and increases to a large enough value that overwhelms the inertia and viscous effects inside the solid phase. As a result, around the gas-solid and liquid-solid interfaces, the momentum equation Eq.(10) reduces to Darcy’s law (Volleretal1987). Voller and Prakash (VollerPrakash1987) modified the well-known Carman-Kozeny equation (Carman1997), by adding a small constant, denoted by here, at the denominator to avoid division by zero, and obtained in , where should be a large number and is the volume fraction of the liquid phase of “” in as a reminder. Such a definition of has been popularly used, e.g., in (Brentetal1988; Shmuelietal2010; RoslerBruggemann2011; Panwisawasetal2017; VogelThess2019). It should be noted that Voller and Prakash (VollerPrakash1987) did not consider the gas phase and formulated in terms of the volume fraction of the liquid phase. Therefore, their formulation is only applicable in . To obtain in , we directly use the volume fraction of the solid phase in as the dependent variable of , i.e., , as shown in Eq.(10). Another popular option of enforcing is to assign a large viscosity inside the solid phase, e.g., in (Volleretal1987; Yanetal2017; Yanetal2018). Although it looks simpler, our numerical tests in Section 4.1.5 show that this is not an effective choice and adding a drag force is recommended.

2.1.6 The energy equation

The enthalpy of the multiphase mixture is , and only the heat conductivity is considered in the present study. As a result, we obtain the following energy equation:

(11)

Here, and denote the volumetric heat and the heat conductivity, respectively, and they are computed from Eq.(5). is the heat source, and is neglected unless otherwise specified.

The terms in the bracket on the left-hand side of Eq.(11) represents the effect of the phase change to the enthalpy. Therefore, without the phase change, i.e., and or and , these terms should disappear. To achieve this goal, in Eq.(11) needs to be replaced by or , due to the consistency of volume fraction conservation. As a result, the terms in the bracket in Eq.(11) after the modification are identical to those on the left-hand side of Eq.(4).

Further, the physical energy equation should admit isothermal solutions (or temperature equilibrium) when the phase change does not happen. Plugging in Eq.(11), we obtain on the left-hand side and on the right-hand side. However, is not zero with computed from Eq.(5). Therefore, the isothermal solution, i.e., , is not admissible by Eq.(11). Following the derivation from the consistency of mass conservation in Section 2.1.4 and replacing with , we can show that is zero, where is the flux of volumetric heat defined in Eq.(5) as well, by noticing that the velocity is divergence-free when there is no phase change. Therefore, on the left-hand side of Eq.(11), needs to be replaced by . Combining the above modifications to Eq.(11), we obtain the consistent energy equation:

(12)

2.2 Properties

Eq.(1), Eq.(4), Eq.(5), Eq.(8), Eq.(10), and Eq.(12) complete the consistent and conservative model for thermo-gas-liquid-solid flows including liquid-solid phase change, and the proposed model honors many physical properties. From Eq.(9), Eq.(10), and Eq.(12), it is obvious that the mass and enthalpy of the multiphase mixture are conserved, and the momentum (neglecting the gravity) is conserved without the appearance of the solid phase, by noticing that is equivalent to , see (Jametetal2002; Jacqmin1999; Shen2011; Huangetal2020N), and , in this circumstance. When there is the solid phase, the momentum is not necessarily conserved due to the no-slip boundary condition at the solid boundary. The Galilean invariance is also satisfied by the proposed model. The proof is straightforward, using the Galilean transformation, and examples are available in (Huangetal2020; Huangetal2020NPMC). The subtle part of the proof is related to the left-hand side of the momentum equation, and we need to emphasize that the consistency conditions are playing a critical role there.

More importantly, the proposed model is reduction consistent with (i) the isothermal consistent and conservative Phase-Field method for two-phase incompressible flows in (Huangetal2020) when the solid phase is absent and the initial homogeneous temperature is larger than the melting temperature, (ii) the isothermal fictitious domain Brinkman penalization (FD/BP) method for fluid-structure interaction in (Angotetal1999; BergmannIollo2011) when the liquid phase is absent and the initial homogeneous temperature is lower than the melting temperature, and (iii) the Phase-Field model of solidification in (Boettingeretal2002) when the gas phase and flow are absent and the material properties of the liquid and solid phases are matched (except the thermal conductivity).

Theorem 2.1.

The proposed model in Section 2.1 is consistent with the isothermal consistent and conservative Phase-Field method for two-phase incompressible flows in (Huangetal2020).

Proof.

Given and at , as already analyzed in Section 2.1.2 and Section 2.1.6, we obtain and therefore and from Eq.(4), and from Eq.(12), at . As a result, the velocity is divergence-free from Eq.(8), and the last term on the right-hand side of the Cahn-Hilliard equation Eq.(1) vanishes. From Eq.(5), we obtain and , showing that the contribution of disappears. Finally in Eq.(10), becomes , and due to . Therefore, the temperature remains its homogeneous initial value, and the simplified system from the proposed model in Section 2.1 with the given condition is equivalent to the consistent and conservative Phase-Field method for two-phase incompressible flows in (Huangetal2020), by noticing that is the order parameter of the Cahn-Hilliard equation in (Huangetal2020). ∎

Theorem 2.2.

The proposed model in Section 2.1 is consistent with the isothermal fictitious domain Brinkman penalization (FD/BP) method for fluid-structure interaction in (Angotetal1999; BergmannIollo2011).

Proof.

Given and at , as already analyzed in Section 2.1.2 and Section 2.1.6, we obtain and therefore and from Eq.(4), and from Eq.(12), at . As a result, the velocity is divergence-free from Eq.(8), and the last term on the right-hand side of the Cahn-Hilliard equation Eq.(1) vanishes. In Eq.(10), becomes , and with . In this case, the gas phase should be understood as an arbitrary incompressible fluid and the solid phase represents the fictitious domain, where the material properties are the same as the fluid ones without loss of generality. As a result, we obtain and from Eq.(5). Therefore, the temperature remains its initial homogeneous value, and the simplified system from the proposed model in Section 2.1 with the given condition is equivalent to the FD/BP method for fluid-structure interaction in (Angotetal1999; BergmannIollo2011). Notice that is defined proportional to in (Angotetal1999; BergmannIollo2011), different from the one in Eq.(10), and the level-set method, instead of the Cahn-Hilliard equation, is used in (BergmannIollo2011) for the volume fraction of the fictitious domain (or the solid phase). ∎

Theorem 2.3.

The proposed model in Section 2.1 is consistent with the Phase-Field model of solidification in (Boettingeretal2002).

Proof.

Given and at , we have and therefore from Eq.(1), which implies and at . Further requiring that the material properties of the liquid and solid phases of “” are identical, we obtain and from Eq.(8) and Eq.(10), respectively. Therefore, we obtain at , and all the convection terms are dropped. Putting all these to the phase change equation Eq.(4) and the energy equation Eq.(12), they become the same as those in (Boettingeretal2002), except that is replaced with in the present study. ∎

Remark: In the proofs of Theorem 2.1, Theorem 2.2, and Theorem 2.3, the given conditions are assumed to be true at

in the entire domain for convenience. Actually, we only need those conditions to be true locally at any moment, and Theorem

2.1, Theorem 2.2, and Theorem 2.3 will again be valid. In other words, the proposed model in Section 2.1 will automatically reduce to the corresponding multiphase models whenever one of the phases is locally absent.

3 Discretization of the governing equations

The numerical procedure to solve the proposed model in Section 2.1 is introduced in this section. The differential operators are discretized with the conservative finite difference method as those in (Huangetal2020), such that the convection terms are discretized by the 5th-order WENO scheme JiangShu1996, while the divergence, gradient, and Laplacian operators are approximated by the 2nd-order central difference (FerzigerPeric2001). These discrete operators have been carefully validated in various studies (Huangetal2020; Huangetal2020CAC; Huangetal2020N; Huangetal2020B; Huangetal2020NPMC). Following the notations in (Huangetal2020), the discrete operators are denoted by , and time levels are indicated by superscript. We use to denote the discretization of the time derivative , and to denote the approximation of from previous time levels. Here, , , and in the 1st-order case, while they are , , and in the 2nd order case. Unless otherwise specified, we use the 2nd-order scheme. The major concern is reproducing the physical connection among the governing equations, discussed in Section 2.1, at the discrete level, following the consistency conditions.

First, the Cahn-Hilliard equation Eq.(1) is solved with the convex splitting scheme in (DongShen2012; Huangetal2020), along with the consistent and conservative boundedness mapping (Huangetal2020CAC). Then, the fully-discretized Cahn-Hilliard equation is recovered and rearranged to be

(13)

Therefore, we obtain and the discrete Phase-Field flux after solving the Cahn-Hilliard equation.

Second, we proceed to solve the phase change equation Eq.(4). To preserve the consistency of volume fraction conservation on the discrete level, it should be noted that and in Eq.(13) are inputs to solve the phase change equation. The fully-discretized phase change equation is

(14)

where represents the WENO reconstruction, denotes the linear interpolation, and is linearized around from Taylor expansion.

After solving Eq.(13) and Eq.(14), , , and are obtained from Eq.(5), noticing that is applied to Eq.(5) due to the consistency of reduction, see (Huangetal2020N). Then, the temperature is updated from the following fully-discretized energy equation:

(15)

where represents the WENO reconstruction and denotes the linear interpolation. It should be noted that the terms in the bracket in Eq.(15) are identical to the left-hand side of the fully-discretized phase change equation in Eq.(14). This numerical correspondence is consistent with the derivation in Section 2.1.6.

Finally, the momentum equation Eq.(10) is solved, majorly based on the 2nd-order projection scheme on a collocated grid (Huangetal2020), which has been carefully analyzed and successfully applied to two- and multi-phase problems Huangetal2020; Huangetal2020N; Huangetal2020CAC. Again, , , and in the momentum equation are directly computed from Eq.(5) and the surface tension force is computed from its definition in Eq.(10) with the balanced-force method (Huangetal2020; Francoisetal2006), while the drag force is treated implicitly. As mentioned in Section 2.1.5, should be predominant over either the inertial or viscous effect inside the solid phase. In other words, from Eq.(10), should be much larger than (inertia) or (viscous). Here denotes the grid size. To achieve this goal, we set as , and is fixed to be . Therefore, will always be thousand times larger than both the inertial and viscous effects inside the solid phase, regardless of the numerical setup or material properties. The divergence of the velocity at the new time level, which appears in the projection scheme, is determined from Eq.(8), i.e.,

(16)

where is obtained from its definition in Eq.(4). On the continuous level, as discussed in Section 2.1.4, with and from Eq.(5) and in Eq.(8), Eq.(9) is implied. However, this is not necessarily true after discretization when the phase change happens and the densities of the liquid and solid phases are not the same. As a result, the consistency of mass conservation and consistency of mass and momentum transport are violated. In order to remedy this issue, a momentum source is added to the momentum equation, where is the residual of the fully-discretized mass conservation equation, i.e.,

(17)

It should be noted that only appears on the discrete level due to discretization errors, and it is exactly zero away from the liquid-solid interface.

Following the above steps, the physical connections among different parts of the proposed model are correctly captured at the discrete level. With similar analyses to those in the proofs of Theorem 2.1, Theorem 2.2, and Theorem 2.3 in Section 2.2, one can easily show that those theorems remain intact on the discrete level. This will be numerically validated in Section 4.1.

4 Results

In this section, various numerical tests are performed to validate and demonstrate the proposed model in Section 2.1. Then, the predictions from the proposed model are compared to experimental data and other simulations. Finally, two challenging setups are performed to illustrate the capability of the proposed model. The formal order of accuracy and the conservation property of the discrete operators in Section 3 have been carefully validated in (Huangetal2020; Huangetal2020CAC; Huangetal2020N) and therefore not repeated here. Unless otherwise specified, the initial velocity is zero, and and are set, where denotes the grid/cell size.

4.1 Validation

We first validate the theorems in Section 2.2 with problems having analytical solutions. Then, the correspondence of the mass conservation and the volume change of the phase change material “” is illustrated. Finally, the effectiveness of the proposed surface tension and drag forces is demonstrated.

4.1.1 Large-Density-Ratio advection

We consider a large-density-ratio advection problem to validate Theorem 2.1 where the solid phase is absent and the temperature is above the melting temperature. The unit domain considered is doubly periodic. A circular drop, whose density is , is initially at the center of the domain with a radius , surrounded by a gas whose density is . The specific heats of the phases are , and the material properties of the liquid phase are shared with the solid phase. The viscosity, heat conduction, surface tension, and gravity are neglected. Other parameters are , , , , and . The solid phase is initially absent, i.e., , and therefore we have initial being . The initial velocity and temperature are and , respectively. The domain is discretized by grid cells, and the time step is determined from .

From Theorem 2.1, the solid phase remains absent, i.e., or at , the temperature remains homogeneous, i.e., at , and the two-phase flow solution is produced, with the above setup. Expected results are obtained and shown in Fig.2. In this setup, the circular drop is translated by the homogeneous velocity. Therefore, there should not be any changes to the shape of the drop and the velocity. From Fig.2 a), we observe that the drop correctly returns to its initial location at , without any deformation, and the streamlines at remain straight and homogeneous. Quantitatively, the difference of the velocity from its initial value is the round-off error, as shown in Fig.2 b). Fig.2 b) also shows that Theorem 2.1 is true due to and at .

It is worth mentioning that the density ratio in this case is , and there are no physical effects, e.g., viscosity and thermal conduction, to homogenize the solution. Without satisfying the consistency conditions, the drop will suffer from unphysical deformations, and the velocity and order parameter will become fluctuating, which are observed in (Huangetal2020; Huangetal2020CAC; Huangetal2020N; Huangetal2020NPMC). The same will happen to the temperature if Eq.(11) is applied, instead of the proposed Eq.(12) that satisfies the consistency conditions, as analyzed in Section 2.1.6.

Figure 2: Results of the large-density-ratio advection. a) Streamlines at and interface of the drop () at and . Blue arrow lines: streamlines at , Black solid line: interface at , Red dashed line: interface at . The drop returns to its original location at without deformation. b) norms of , , , and versus time. The velocity preserves its initial value, the solid phase remains absent, and the temperate equilibrium is maintained exactly, which validate Theorem 2.1.

4.1.2 Couette flow

A Couette flow problem is solved to validate Theorem 2.2 where the liquid phase is absent and the temperature is below the melting temperature. The unit domain considered is periodic along the axis while is no-slip along the axis. Both the top and bottom boundaries are adiabatic but the top one is moving with a unit horizontal velocity, i.e., . The solid phase is at the bottom below , while the gas phase fills the rest of the domain. The input parameters are , , , , , , , , , and . The liquid phase is initially absent, i.e., , and therefore we have . The initial temperature is which is lower than the melting temperature . The domain is discretized by grid cells, and the time step is determined by .

The above setup is equivalent to the following Couette flow:

(18)

where , , , and . The exact solution of Eq.(18) is

(19)

derived from separation of variables. Theorem 2.2 implies that the liquid phase remains absent, i.e., or at , the temperature remains homogeneous, i.e., at , and the solution in Eq.(19) is produced (or approximated) by the FD/BP FSI formulation, with the present setup. Expected results are obtained and shown in Fig.3. First in Fig.3 a), the profile of the solid fraction at overlaps the one at , representing that the solid phase is stationary. Moreover, the profiles of the solution from the proposed model are indistinguishable from the exact solution in Eq.(19), noticing that the summation in Eq.(19) is from to . Theorem 2.2 is true, as shown in Fig.3 b) that and at . In addition, the unidirectional condition, which is required to obtain Eq.(18), is also demonstrated, as shown in Fig.3 b) that at .

Figure 3: Results of the Couette flow. a) Profile of at and , and profiles of and at , , , , and . Here, is the exact solution of the Couette flow in Eq.(19), and the numerical predictions and analytical solutions are overlapped. b) norms of , , and versus time. The liquid phase remains absent, the temperature equilibrium is preserved, and the unidirectional flow condition is produced, which validate Theorem 2.2.

4.1.3 Stefan problem

The Stefan problem is performed to validate Theorem 2.3 where both the gas phase and the flow are absent and the liquid and solid phases have matched material properties (except the thermal conductivities). The setup in (Javierreetal2006) is followed. The unit domain is periodic along the direction. Both the top and bottom boundaries are free-slip and adiabatic. The material properties are: , , , , , , , and . Other parameters are , , and , where , , and , the same as those in (Javierreetal2006). Initially, the liquid-solid interface is at , above which there is the solid phase having a temperature , while below which there is the liquid phase having a temperature . The domain is discretized by grid cells, and the time step is .

The above setup is to produce the following Stefan problem (Javierreetal2006) (based on ):

(20)

where is the location of the liquid-solid interface. Eq.(20) has an analytical self-similar solution (James1987; Javierreetal2006):

(21)

Numerical results are compared to the exact solution Eq.(21), and shown in Fig.4. The interface location is specified as the contour of from the numerical results. Both the temperature and interface location agree with the exact solution very well. Minor discrepancy is observed near the domain boundary at in Fig.4 a) because in practice the domain is not infinite. Moreover, the norms of , , and are on the orders of , , and , respectively, at the end of the simulation, which demonstrates Theorem