1 Introduction
LiquidSolid 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 highfidelity models to further understand the complex physics and dynamics, and accurately predict the behavior of materials in order to produce highquality products. Such a problem includes many challenging factors, such as a wide range of material properties, evolution of the liquidsolid interface due to solidification/melting, deformation of the gasliquid interface due to fluid (gas/liquid) motions and surface tension, heat transfer that drives the phase change, and fluidstructure interaction between the fluid and solid, and all these factors are coupled and can be influential. In the present study, we call the problem thermogasliquidsolid flows with liquidsolid phase change.
In spite of the complexity of the problem, it can be separated into two basic problems which are the twophase incompressible flow and the solidification with convection. Several studies focused on these individual problems. For twophase incompressible flow, the fronttracking method UnverdiTryggvason1992; Tryggvasonetal2001, the levelset method OsherSethian1988; Sussmanetal1994; SethianSmereka2003; Gibouetal2018, the conservative levelset method OlssonKreiss2005; Olssonetal2007; ChiodiDesjardins2017, the volumeoffluid (VOF) method HirtNichols1981; ScardovelliZaleski1999; OwkesDesjardins2017, the THINC method Xiaoetal2005; Iietal2012; XieXiao2017; Qianetal2018, and the PhaseField (or DiffuseInterface) 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 wellbalanced surface tension model AbuAlSaud2018, and the PhaseField method derived from the energy balance or the leastaction principle Jacqmin1999; Yueetal2004 have been developed to model the surface tension, and a balancedforce 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 enthalpyporosity 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 fluidstructure interaction (Angotetal1999; BergmannIollo2011). The most popularly used drag force in the enthalpyporosity technique is modified from the CarmanKozeny equation (Carman1997) for the porous medium. Another popular method for solidification is the PhaseField method (Boettingeretal2002; Chen2002; Echebarriaetal2004; KimKim2005; Tanetal2011; Jietal2018; Luetal2018), where the liquidsolid interface has a small but finite thickness. Different from the enthalpyporosity technique using an algebraic relation, the PhaseField 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 PhaseField method with the hydrodynamics, the melt convection can be modeled (Nestleretal2000; Beckermannetal1999; ChenYang2019; ZhangYang2020). Other methods for modeling the liquidsolid phase change are reviewed in (SalcudeanAbdullah1988; Samarskiietal1993; Voller1996; HuArgyropoulos1996; Dutiletal2011; DhaidanKhodadadi2015; Sultanaetal2018).
Most existing models for the thermogasliquidsolid flows with liquidsolid phase change follow a similar procedure: The deforming gasliquid interface is located by a certain interface capturing method and the continuous surface tension force is added, while the enthalpyporosity technique is directly applied without any further changes to adapt to the appearance of a new phase. The volumeoffluid method is the most popular choice and is used, e.g., in (Shmuelietal2010; Saldi2012; Kimetal2013; Yanetal2017; Panwisawasetal2017; VogelThess2019; Heetal2020). The levelset and conservative level set methods are recently used in (Yanetal2018) and (Linetal2020), respectively. Some additional physics are introduced to the models, e.g., the thermocapillary effect (Panwisawasetal2017; Yanetal2018; Heetal2020; Linetal2020) and recoil pressure (Panwisawasetal2017; Heetal2020; Linetal2020). Another recent model (Zhangetal2020) follows the same strategy but instead uses the PhaseField model in (Ramirezetal2004) for anisotropic solidification and the conservative PhaseField method (ChiuLin2011) as the interface capturing method. In spite of its widespread applications, such a wellaccepted modeling strategy has the following critical issues. (i) The volume fractions of the phases are ambiguously defined. In the models applying the enthalpyporosity 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 gasliquid interface will falsely appear at the gassolid 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 bounceback scheme near the gassolid interface is employed since the gassolid interface is unable to be effectively labeled by the model, but details of the bounceback 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 divergencefree, 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 divergencefree without evaporation. As a result, the divergencefree 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 wellknown that the socalled 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 twophase flow. In the present study, those critical issues are properly addressed and a consistent and conservative PhaseField model is developed for the thermogasliquidsolid flows with liquidsolid 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 multifold 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 PhaseField method for incompressible twophase flows (Huangetal2020) when the solid phase is absent, the fictitious domain Brinkman penalization (FD/BP) method for fluidstructure interactions (Angotetal1999; BergmannIollo2011) when the liquidphase is absent, and the PhaseField 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 nonzero, and therefore the volume change during solidification/melting is captured.

The momentum is conserved when the solid phase is absent, and the noslip 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 gasliquidsolid mixture, which greatly improves the robustness of the model for largedensityratio 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/solidstate 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 liquidsolid 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 thermogasliquidsolid flows including liquidsolid 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/solidstate.

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 thermogasliquidsolid problems locally are gasliquid, gassolid, or liquidsolid problems and can be isothermal. Therefore, all these circumstances should be admissible by a thermogasliquidsolid 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 PhaseField 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 CahnHilliard equation
The interfacial dynamics of “” and “” is modeled by the CahnHilliard equation (CahnHilliard1958) with convection:
(1)  
Here, is the order parameter of the CahnHilliard 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 gasliquid interface . is the doublewell potential, having minimums at or , and is the derivative of with respect to . is the PhaseField flux of , including both the convection and diffusion fluxes in the CahnHilliard equation. Unless otherwise specified, the homogeneous Neumann boundary condition is applied.
The CahnHilliard equation Eq.(1) is derived from the GinzburgLandau free energy functional using the gradient flow, and the chemical potential in Eq.(1) is the functional derivative of the GinzburgLandau free energy functional with respect to the order parameter . The CahnHilliard equation has been widely used in modeling twophase 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 AllenCahn PhaseField 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 AllenCahn dynamics, and is considered as the volume fraction of the liquid phase of “” in . is the mobility of . and are the GibbsThomson and linear kinetic coefficients, respectively, in the GibbsThomson equation of the liquidsolid interface. is the mixing energy density of , which is related to the thickness of the liquidsolid interface . is the latent heat of the liquidsolid 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 doublewell 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 GibbsThomson equation, and details are available in (Boettingeretal2002; Beckermannetal1999; AllenCahn1979). The first two terms in the parentheses on the righthand 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 liquidsolid 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 righthand 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 lefthand side and on the righthand side, given and , and these two sides are equal to each other from the CahnHilliard 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 PhaseField flux of .
2.1.3 The volume fractions and material properties
Based on the CahnHilliard 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.
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 CahnHilliard 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 liquidsolid phase change. As a result, the mass conservation equation of the model is
(9) 
Therefore, the mass of the gasliquidsolid 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 gasliquid, gassolid, and liquidsolid 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 gasliquid interface, while is the drag force, modeling the noslip boundary condition at the gassolid and liquidsolid 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 righthand side of Eq.(10)) and the Galilean invariance, as analyzed in (Huangetal2020).
To model the surface tension at the gasliquid interface, we apply the commonly used PhaseField formulation , which can be derived from either the energy balance Jacqmin1999; Huangetal2020CAC or the leastaction principle Yueetal2004; Shen2011. However, is activated at all “” interfaces, including both the gasliquid and gassolid interfaces. To remove its contribution at the gassolid interface, is multiplied by so that the surface tension force remains to be at the gasliquid 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 thermocapillary effect has not been considered in the present study and the surface tension is treated as a constant.
To model the noslip boundary condition at both the gassolid and liquidsolid 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 gassolid and liquidsolid interfaces, the momentum equation Eq.(10) reduces to Darcy’s law (Volleretal1987). Voller and Prakash (VollerPrakash1987) modified the wellknown CarmanKozeny 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 lefthand 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 lefthand 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 lefthand side and on the righthand 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 divergencefree when there is no phase change. Therefore, on the lefthand 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 thermogasliquidsolid flows including liquidsolid 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 noslip 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 lefthand 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 PhaseField method for twophase 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 fluidstructure interaction in (Angotetal1999; BergmannIollo2011) when the liquid phase is absent and the initial homogeneous temperature is lower than the melting temperature, and (iii) the PhaseField 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 PhaseField method for twophase 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 divergencefree from Eq.(8), and the last term on the righthand side of the CahnHilliard 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 PhaseField method for twophase incompressible flows in (Huangetal2020), by noticing that is the order parameter of the CahnHilliard 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 fluidstructure 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 divergencefree from Eq.(8), and the last term on the righthand side of the CahnHilliard 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 fluidstructure interaction in (Angotetal1999; BergmannIollo2011). Notice that is defined proportional to in (Angotetal1999; BergmannIollo2011), different from the one in Eq.(10), and the levelset method, instead of the CahnHilliard 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 PhaseField 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 5thorder WENO scheme JiangShu1996, while the divergence, gradient, and Laplacian operators are approximated by the 2ndorder 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 1storder case, while they are , , and in the 2nd order case. Unless otherwise specified, we use the 2ndorder 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 CahnHilliard equation Eq.(1) is solved with the convex splitting scheme in (DongShen2012; Huangetal2020), along with the consistent and conservative boundedness mapping (Huangetal2020CAC). Then, the fullydiscretized CahnHilliard equation is recovered and rearranged to be
(13) 
Therefore, we obtain and the discrete PhaseField flux after solving the CahnHilliard 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 fullydiscretized 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 fullydiscretized 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 lefthand side of the fullydiscretized 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 2ndorder projection scheme on a collocated grid (Huangetal2020), which has been carefully analyzed and successfully applied to two and multiphase 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 balancedforce 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 fullydiscretized 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 liquidsolid 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 LargeDensityRatio advection
We consider a largedensityratio 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 twophase 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 roundoff 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.
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 noslip 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 .
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 freeslip and adiabatic. The material properties are: , , , , , , , and . Other parameters are , , and , where , , and , the same as those in (Javierreetal2006). Initially, the liquidsolid 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 liquidsolid interface. Eq.(20) has an analytical selfsimilar 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
Comments
There are no comments yet.