A consistent and conservative Phase-Field model for thermo-gas-liquid-solid flows including liquid-solid phase change
AA consistent and conservative Phase-Field model forthermo-gas-liquid-solid flows including liquid-solid phase change
Ziyang Huang ∗ , Guang Lin † , and Arezoo M. Ardekani ‡ School of Mechanical Engineering, Purdue University, West Lafayette, IN 47907, USA Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA (February 13, 2021)
Abstract
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 theconsistency conditions and exactly reduces to the consistent and conservative Phase-Field method forincompressible 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 massconservation, defines the volume fractions of individual phases unambiguously, and therefore captures thevolume change due to phase change. The momentum is conserved when the solid phase is absent, but itchanges when the solid phase appears due to the no-slip condition at the solid boundary. The proposedmodel also conserves the energy, preserves the temperature equilibrium, and is Galilean invariant. Anovel continuous surface tension force to confine its contribution at the gas-liquid interface and a dragforce modified from the Carman-Kozeny equation to reduce solid velocity to zero are proposed. The issueof initiating phase change in the original Phase-Field model of solidification is addressed by physicallymodifying the interpolation function. The corresponding consistent scheme is developed to solve themodel, and the numerical results agree well with the analytical solutions and the existing experimentaland numerical data. Two challenging problems having a wide range of material properties and complexdynamics are conducted to demonstrate the capability of the proposed model.
Keywords:
Consistent model; Phase change; Solidification/Melting; Multiphase flow; Fluid-Structureinteraction; Phase-Field;
Liquid-Solid phase change (or solidification/melting) and its interaction with the surrounding air are ubiq-uitous in various natural and/or industrial processes, e.g., latent thermal energy storage (LTES) systems[80, 88], welding [65, 13, 101, 74], casting [19, 36], and additive manufacturing (AM) [64, 32, 53]. This mo-tivates researchers to develop physical and high-fidelity models to further understand the complex physicsand 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 ofthe 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 in-teraction between the fluid and solid, and all these factors are coupled and can be influential. In the presentstudy, 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 ∗ Email: [email protected] † Email: [email protected] ; Corresponding author ‡ Email: [email protected] ; Corresponding author a r X i v : . [ phy s i c s . c o m p - ph ] F e b roblems. For two-phase incompressible flow, the front-tracking method [87, 85], the level-set method[62, 82, 77, 30], the conservative level-set method [60, 61, 17], the volume-of-fluid (VOF) method [33, 76, 63],the THINC method [93, 42, 95, 68], and the Phase-Field (or Diffuse-Interface) method [3, 43, 79, 41] havebeen developed to locate different phases. The smoothed surface stress method [31], the continuous surfaceforce (CSF) [8], the ghost fluid method (GFM) [25, 51], the conservative and well-balanced surface tensionmodel [57], and the Phase-Field method derived from the energy balance or the least-action principle [43, 98]have been developed to model the surface tension, and a balanced-force method [28] is proposed to incorporatethe surface tension model to the fluid motion. Interested readers should refer to [67, 86, 79, 56, 66]. For thesolidification problem with convection, one of the most popular methods is the enthalpy-porosity technique[89, 90, 9, 91, 71], which can be implemented in a fixed grid. The liquid fraction of the phase change materialis algebraically related to the local temperature. Therefore, the enthalpy change due to phase change can beevaluated and becomes a source in the energy equation. A drag force proportional to the velocity is added tothe momentum equation to stop the solid motion, which is the same idea as the fictitious domain Brinkmanpenalization (FD/BP) method for the fluid-structure interaction [4, 6]. The most popularly used drag forcein the enthalpy-porosity technique is modified from the Carman-Kozeny equation [12] for the porous medium.Another popular method for solidification is the Phase-Field method [7, 15, 24, 48, 83, 46, 54], where theliquid-solid interface has a small but finite thickness. Different from the enthalpy-porosity technique usingan algebraic relation, the Phase-Field method introduces an additional equation to govern the evolutionof the liquid fraction of the phase change material, and therefore is flexible to include more complicatedphysics, e.g., anisotropy and solute transport in an alloy. After coupling the Phase-Field method with thehydrodynamics, the melt convection can be modeled [59, 5, 14, 100]. Other methods for modeling theliquid-solid phase change are reviewed in [73, 75, 92, 35, 23, 20, 81].Most existing models for the thermo-gas-liquid-solid flows with liquid-solid phase change follow a similarprocedure: The deforming gas-liquid interface is located by a certain interface capturing method and thecontinuous surface tension force is added, while the enthalpy-porosity technique is directly applied withoutany further changes to adapt to the appearance of a new phase. The volume-of-fluid method is the mostpopular choice and is used, e.g., in [80, 74, 49, 97, 64, 88, 32]. The level-set and conservative level set methodsare recently used in [96] and [53], respectively. Some additional physics are introduced to the models, e.g., thethermo-capillary effect [64, 96, 32, 53] and recoil pressure [64, 32, 53]. Another recent model [99] follows thesame strategy but instead uses the Phase-Field model in [70] for anisotropic solidification and the conservativePhase-Field method [18] as the interface capturing method. In spite of its widespread applications, sucha well-accepted modeling strategy has the following critical issues. (i) The volume fractions of the phasesare ambiguously defined. In the models applying the enthalpy-porosity technique, the liquid fraction ismeaningful only inside the phase change material, while it directly appears in the energy equation definedin the entire domain including the gas phase. As a result, the meaningless value of the liquid fraction in thegas phase is also counted in the energy equation. Another example is in [99], where two liquid fractions aredefined for the same liquid phase, one from the solidification model and the other from the interface capturingmethod. 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 forcescan appear at the wrong location because the volume fractions, which are not clearly defined, are neededto compute the forces. Based on the formulations, e.g., in [74, 64, 96, 53, 32], the surface tension at thegas-liquid interface will falsely appear at the gas-solid interface, and the drag force will falsely appear in thegas phase when the local temperature is lower than the solidus temperature. Some artificial operations needto be added but they have seldom been mentioned in the literature. An exception is in [99] where a bounce-back scheme near the gas-solid interface is employed since the gas-solid interface is unable to be effectivelylabeled by the model, but details of the bounce-back scheme are not provided. (iii) Physical principles canbe violated, depending on the material properties.
The most obvious example is the mass conservation. In,e.g., [71, 64, 97, 96, 99, 32, 49], the velocity is divergence-free, implying that both the volumes of the gasand phase change material will not change. This restricts the application of the model to problems havingmatched liquid and solid densities. However, the problems studied in [49, 53] are outside that category, andtherefore 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 [53] is from evaporation, not solidification/melting, and thevelocity is divergence-free without evaporation. As a result, the divergence-free velocity is contradicting themass conservation. Although the studies in [78, 80, 34, 88] captures the volume change, detail formulations,2.e., the divergence of the velocity, are not provided. Other physical principles, e.g., the momentum andenergy conservation and the Galilean invariance of the models, have never been examined. (iv) The problemof interest often has a large density ratio , which can be O (10 ) between the gas and liquid. It has beenwell-known that the so-called consistent method [72, 10, 16, 63, 69, 58, 94, 41, 39] needs to be implementedto produce physical results for multiphase flows, while this has never been considered in the existing modelsfor 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 presentstudy, those critical issues are properly addressed and a consistent and conservative Phase-Field model isdeveloped for the thermo-gas-liquid-solid flows with liquid-solid phase change. All the dependent variablesare 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 physicalproperties: • In deriving the proposed model, several consistency conditions proposed in our previous works [41, 39,38, 40, 37] are considered. Our earlier works do not include phase change or temperature variationsin the fluids, and the present study is the first implementation of the consistency conditions to phasechange problems, which further demonstrates their generality. • The proposed model exactly recovers the consistent and conservative Phase-Field method for incom-pressible two-phase flows [41] when the solid phase is absent, the fictitious domain Brinkman penaliza-tion (FD/BP) method for fluid-structure interactions [4, 6] when the liquid-phase is absent, and thePhase-Field model of solidification of a pure material [7] when the gas phase is absent. • The proposed model defines the volume fractions of individual phases unambiguously and ensures theirsummation to be unity everywhere. The local mass conservation is strictly satisfied, from which thedivergence of the velocity is non-zero, and therefore the volume change during solidification/melting iscaptured. • The momentum is conserved when the solid phase is absent, and the no-slip condition at the solidboundary results in the momentum change. The energy conservation and Galilean invariance are alsosatisfied by the proposed model. • The momentum transport is consistent with the mass transport of the gas-liquid-solid mixture, whichgreatly improves the robustness of the model for large-density-ratio problems and avoids unrealisticinterface deformation. • Isothermal (or temperature equilibrium) solutions are admissible, thanks to satisfying the consistencyconditions, which prevents producing any fictitious fluctuations of the temperature. • Novel continuous surface tension and drag force models are proposed, which are activated only atproper locations. The interpolation function in the solidification model [7] is modified to include thecapability of initiating phase change when there is only the liquid/solid-state of the phase changematerial.These physical properties of the proposed model are independent of the material properties.The proposedmodel is validated and its capability is demonstrated using a consistent numerical scheme. The proposedmodel includes all the basic ingredients and challenging aspects of the problem, and additional physics canbe incorporated conveniently following the same framework.The rest of the paper is organized as follows. In Section 2, the proposed model and its properties areelaborated 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, andchallenging problems are simulated to demonstrate the capability of the proposed model. In Section 5, thepresent study is concluded. 3
Governing equations
The problem considered includes two materials, which are a gas “ G ” and a phase change material “ M ”experiencing solidification or melting. Therefore, in the entire domain Ω, there are three phases: the gasphase including only “ G ”, and the liquid and solid phases of “ M ”. The liquid-solid phase change is drivenby temperature. The part of Ω occupied by “ G ” is denoted by Ω G , and similarly, we use Ω M , Ω LM , and Ω SM to denote the domains occupied by “ M ”, the liquid phase, and the solid phase of “ M ”, respectively. As aresult, we have Ω = Ω G ∪ Ω M = Ω G ∪ Ω LM ∪ Ω SM . In addition, boundaries of the domains are denoted with“ ∂ ” in front of the corresponding domains, for example, ∂ Ω is the boundary of Ω. The material propertiesof the gas phase “ G ” and the liquid and solid phases of “ M ” are assumed to be constant and denoted by β G , β LM , and β SM , respectively, where β can be density ρ , viscosity µ , specific heat C p , and heat conductivity κ .The proposed consistent and conservative model for thermo-gas-liquid-solid flows including liquid-solidphase change is elaborated in Section 2.1. Then in Section 2.2, the physical properties and the relations ofthe proposed model to some other multiphase models are analyzed. When the proposed model is derivedin Section 2.1, several consistency conditions will be applied so that the proposed model is able to producephysical results. The definitions of the consistency conditions are: • Consistency of reduction:
The multiphase system should be able to recover the corresponding systemsincluding fewer phases. • Consistency of volume fraction conservation:
The phase change equation should be consistent withthe governing equation for the volume fraction of “ M ”, when “ M ” is in a fully liquid/solid-state. • Consistency of mass conservation:
The mass conservation equation should be consistent with thegoverning equation for the volume fraction of “ M ”, the phase change equation, and the density ofthe multiphase mixture. The mass flux and the divergence of the velocity in the mass conservationequation should lead to a zero mass source. • Consistency of mass and momentum transport:
The momentum flux in the momentum equation shouldbe computed as a tensor product between the mass flux and the velocity, where the mass flux shouldbe identical to the one in the mass conservation equation.These consistency conditions are generalized from their correspondences for isothermal multimaterial incom-pressible flows [41, 39, 38, 40, 37] 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 consistencyof reduction . To physically understand the meaning of the rest of the consistency conditions, one needs tofirst realize that the Phase-Field method includes relative motions of the materials, modeled as a diffusiveprocess. The consistency conditions are general principles to incorporate the mass, momentum, and energytransports due to the relative motions of the materials into the conservation equations.
The interfacial dynamics of “ G ” and “ M ” is modeled by the Cahn-Hilliard equation [11] with convection: ∂ϕ∂t + ∇ · ( u ϕ ) = ∇ · ( M ϕ ∇ ξ ϕ ) + ϕ ∇ · u in Ω , (1) ξ ϕ = λ ϕ (cid:18) η ϕ g (cid:48) ( ϕ ) − ∇ ϕ (cid:19) ,λ ϕ = 3 √ η ϕ σ, g ( ϕ ) = ϕ (1 − ϕ ) , m ϕ = u ϕ − M ϕ ∇ ξ ϕ , n · ∇ ξ ϕ = n · ∇ ϕ = 0 at ∂ Ω . ϕ is the order parameter of the Cahn-Hilliard equation and is considered as the volume fraction of “ M ”in Ω. u is the velocity, whose divergence can be nonzero. M ϕ and ξ ϕ are the mobility and chemical potentialof ϕ , respectively. λ ϕ is the mixing energy density of ϕ , which is proportional to both the thickness of “ G - M ”interface η ϕ and the surface tension at the gas-liquid interface σ . g ( ϕ ) is the double-well potential, havingminimums at ϕ = 0 or ϕ = 1, and g (cid:48) ( ϕ ) is the derivative of g ( ϕ ) with respect to ϕ . m ϕ is the Phase-Fieldflux of ϕ , including both the convection and diffusion fluxes in the Cahn-Hilliard equation. Unless otherwisespecified, the homogeneous Neumann boundary condition is applied.The Cahn-Hilliard equation Eq.(1) is derived from the Ginzburg-Landau free energy functional using the H − 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 beenwidely used in modeling two-phase incompressible flows, e.g., in [43, 21, 1, 41], and its derivation has alreadybeen given in many references, e.g., in [79, 26, 98]. In the present study, we consider the Allen-Cahn Phase-Field model for the solidification of a pure material[7] and the convection term is added: ∂φ∂t + ∇ · ( u φ ) = − M φ (cid:32) λ φ η φ g (cid:48) ( φ ) − λ φ ∇ φ + ρ LM LT M p (cid:48) ( φ )( T M − T ) (cid:33) + φ ∇ · u in Ω M , (2) λ φ = 3 √ ρ LM LT M Γ φ η φ , M φ = µ φ Γ φ λ φ ,p ( φ ) = φ (6 φ − φ + 10) , n · ∇ φ = 0 at ∂ Ω M . Notice that Eq.(2) is defined in Ω M . Here, φ is the order parameter of the Allen-Cahn dynamics, and isconsidered as the volume fraction of the liquid phase of “ M ” in Ω M . M φ is the mobility of φ . Γ φ and µ φ are the Gibbs-Thomson and linear kinetic coefficients, respectively, in the Gibbs-Thomson equation ofthe liquid-solid interface. λ φ is the mixing energy density of φ , which is related to the thickness of theliquid-solid interface η φ . L is the latent heat of the liquid-solid phase change. T is the temperature and T M is the melting temperature of the phase change. p ( φ ) is the interpolation function, monotonically increasingfrom 0 to 1 and having extreme points at φ = 0 and φ = 1, and p (cid:48) ( φ ) is the derivative of p ( φ ) with respectto φ . g ( φ ) is the double-well potential function defined identically to the one in Eq.(1). Unless otherwisespecified, 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 [7, 15, 46, 48, 14, 100]. The model can be derived either thermodynamicallyfrom a free energy functional using the L gradient flow or geometrically from the Gibbs-Thomson equation,and details are available in [7, 5, 2]. 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 witheach other to maintain the thickness of the liquid-solid interface. The effect of the temperature is modeledby the last term in the parentheses.Solving Eq.(2) is very challenging because it is defined in Ω M which is evolving with time. It would bemore convenient if we can obtain an equivalent equation to Eq.(2) but defined in Ω. Therefore, the diffusedomain approach [52] is applied, and we use ϕ , the volume fraction of “ M ”, as an approximation to theindicator function of Ω M whose value is 1 in Ω M but 0 elsewhere. The equivalence of Eq.(2) in Ω is ∂ ( ϕφ ) ∂t + ∇ · ( u ϕφ ) = − M φ (cid:32) λ φ η φ ϕg (cid:48) ( φ ) − λ φ ∇ · ( ϕ ∇ φ ) + ρ LM LT M ϕp (cid:48) ( φ )( T M − T ) (cid:33) + ϕφ ∇ · u in Ω , (3) n · ∇ φ = 0 at ∂ Ω . We directly apply the formulations of the diffuse domain approach in [52] to obtain Eq.(3) from Eq.(2), anddetails of the approach are available in [52]. Two modifications will be performed to Eq.(3) to address thefollowing two issues. 5he first issue is about initiating the phase change. It should be noted that the terms in the parentheseson the right-hand side of Eq.(2), which are the driving forces for the phase change, are nonzero only at0 < φ <
1. In other words, given “ M ” 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 p ( φ ) in Eq.(2) whose extreme points are φ = 0 and φ = 1, i.e., p (cid:48) (0) = p (cid:48) (1) = 0. These extreme points arethe same as the equilibrium states of φ . Defining p ( φ ) in this way, as mentioned in [7], is an improvementfrom using p ( φ ) = φ , e.g., in [45], in the sense that the equilibrium state of φ is always 0 or 1, independentof the temperature. However, defining p ( φ ) = φ preserves the driving force from the temperature when φ = 0 or φ = 1. The equilibrium states of φ should depend on the temperature. If the temperature is largerthan the melting point, the equilibrium state should be φ = 1 (liquid), while it should be φ = 0 (solid) if T < T M . To achieve this property, we propose ˜ p (cid:48) ( φ ) defined in Eq.(4) in the present study, which combinesthe advantages of p ( φ ) = φ (6 φ − φ + 10) in [7] and p ( φ ) = φ in [45], but avoids their disadvantages.Only when φ = 1 (liquid) and T < T M (undercool) or when φ = 0 (solid) and T > T M (overheat), ˜ p (cid:48) ( φ ) is1, so that the effect of the temperature on the phase change is included. In other cases, ˜ p (cid:48) ( φ ) is the same as p (cid:48) ( φ ) with p ( φ ) = φ (6 φ − φ + 10). As a result, the equilibrium state of φ when T < T M is φ = 0 (solid),and it is φ = 1 (liquid) when T > T M .The second issue is about the existence of fully liquid/solid state of “ M ”. For example, given, φ = 1and T > T M , we obtain ∂φ/∂t = 0 from Eq.(2), which implies φ ≡
1. In other words, Eq.(2) admits theexistence of fully liquid state of “ M ”, when the temperature is larger than the melting temperature andthere is no solid phase at the beginning. This property should be inherited by Eq.(3), while this is not thecase. Using the same condition, we are unable to obtain φ ≡ consistency of volume fraction conservation . After comparing Eq.(3) along with φ = 1 and T > T M toEq.(1), we discover that the convection velocity u ϕ in Eq.(3) needs to be replaced with m ϕ . As a result, wehave ∂ϕ/∂t + ∇ · m ϕ on the left-hand side and ϕ ∇ · u on the right-hand side, given φ = 1 and T > T M , andthese two sides are equal to each other from the Cahn-Hilliard equation Eq.(1). Therefore, the existence offully liquid state of “ M ” is admitted by Eq.(3) after the modification. The same is also true for the solidstate of “ M ”.After applying the above mentioned modifications to Eq.(3) to address these two issues, we obtain thephase change equation: ∂ ( ϕφ ) ∂t + ∇ · ( m ϕ φ ) = − M φ ξ φ + ϕφ ∇ · u in Ω , (4) ξ φ = − λ φ ∇ · ( ϕ ∇ φ ) + λ φ η φ ϕg (cid:48) ( φ ) + ρ LM LT M ϕ ˜ p (cid:48) ( φ )( T M − T ) ,λ φ = 3 √ ρ LM LT M Γ φ η φ , M φ = µ φ Γ φ λ φ ,p ( φ ) = φ (6 φ − φ + 10) , ˜ p (cid:48) ( φ ) = , φ = 0 and T (cid:62) T M , , φ = 1 and T (cid:54) T M ,p (cid:48) ( φ ) , else , m ϕφ = m ϕ φ, n · ∇ φ = 0 , at ∂ Ω . Here, ξ φ is called the chemical potential of φ , ˜ p (cid:48) ( φ ) is modified from p (cid:48) ( φ ) to initiate the phase change, and m ϕφ is the Phase-Field flux of ( ϕφ ). Based on the Cahn-Hilliard equation Eq.(1) and the phase change equation Eq.(4), the volume fractions ofthe gas, liquid, and solid phases in the proposed model are defined unambiguously. The volume fraction ofthe gas phase (or “ G ”) in Ω is α G = (1 − ϕ ), and they are α L = ( ϕφ ) and α S = ( ϕ − ϕφ ) for the liquidand solid phases, respectively. The volume fraction of “ M ” in Ω is α M = α L + α S = ϕ . It is clear that α G + α M = α G + α L + α S = 1 is always true in the entire domain.With the volume fractions of the phases in hand, the material properties of the gas-liquid-solid mixture6nd their fluxes are computed as β = β G + ( β SM − β G ) ϕ + ( β LM − β SM )( ϕφ ) , (5) m β = β G u + ( β SM − β G ) m ϕ + ( β LM − β SM ) m ϕφ . Here, β represents a certain material property, e.g., density ρ , and m ϕ and m ϕφ are the Phase-Field fluxesdefined in Eq.(1) and Eq.(4), respectively. To determine the mass transport of the model, the consistency of mass conservation is applied. First, thedensity of the multiphase mixture is obtained following Eq.(5) ρ = ρ G + ( ρ SM − ρ G ) ϕ + ( ρ LM − ρ SM )( ϕφ ) . (6)After combining Eq.(6) with the Cahn-Hilliard equation, Eq.(1), and the phase change equation, Eq.(4), themass of the multiphase mixture is governed by ∂ρ∂t + ∇ · m ρ = ρ ∇ · u − M φ ( ρ LM − ρ SM ) ξ φ , (7)where m ρ is the consistent mass flux defined in Eq.(5). Finally, to obtain a zero mass source in Eq.(7), thedivergence of the velocity should satisfy ∇ · u = M φ ( ρ LM − ρ SM ) ρ ξ φ . (8)Eq.(8) illustrates the volume change due to the liquid-solid phase change. As a result, the mass conservationequation of the model is ∂ρ∂t + ∇ · m ρ = 0 . (9)Therefore, the mass of the gas-liquid-solid mixture is locally conserved, even though the phase changehappens. This is achieved at the expense of changing the volume of “ M ”, as indicated in Eq.(8). On theother hand, if the densities of the liquid and solid phases of “ M ” are the same, the volume of “ M ” remainsthe same, and therefore the divergence of the velocity is zero, which can also be seen in Eq.(8). It should benoted that the mass conservation equation Eq.(9) is not an independent equation in the model. Instead, itis derived from Eq.(5) and Eq.(8), as elaborated in this section. Therefore, we don’t need to explicitly solvethe mass conservation equation Eq.(9). The motion of the phases (or materials) is governed by the momentum equation: ∂ ( ρ u ) ∂t + ∇ · ( m ρ ⊗ u ) = −∇ P + ∇ · [ µ ( ∇ u + ∇ u T )] + ρ g + f s + f d , (10) f s = φξ ϕ ∇ ϕ, f d = A d ( u S − u ) , A d = C d α S (1 − α S ) + e d . Here, P is the pressure, and g is the gravity. As all the gas-liquid, gas-solid, and liquid-solid interfaces areimmersed in Ω, their effects are modeled as volumetric forces, i.e., f s and f d , in the momentum equationEq.(10). f s is the surface tension force, modeling the surface tension at the gas-liquid interface, while f d is the drag force, modeling the no-slip boundary condition at the gas-solid and liquid-solid interfaces byenforcing u = u S in the solid phase. A d is the drag coefficient of f d , and C d and e d are model parameters of A d . u S is the given solid velocity and is set to be zero unless otherwise specified. Note that the momentumis transported with the consistent mass flux m ρ 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 energy7igure 1: Magnitude of the surface tension force | f s | . a) The proposed surface tension force f s = φξ ϕ ∇ ϕ , b)The original surface tension force f s = ξ ϕ ∇ ϕ . The proposed surface tension force shown in a) removes itscontribution at the gas-solid interface, compared to the original formulation shown in b).conservation (when only the pressure gradient is present on the right-hand side of Eq.(10)) and the Galileaninvariance, as analyzed in [41].To model the surface tension at the gas-liquid interface, we apply the commonly used Phase-Field for-mulation ξ ϕ ∇ ϕ , which can be derived from either the energy balance [43, 39] or the least-action principle[98, 79]. However, ξ ϕ ∇ ϕ is activated at all “ G - M ” interfaces, including both the gas-liquid and gas-solidinterfaces. To remove its contribution at the gas-solid interface, ξ ϕ ∇ ϕ is multiplied by φ so that the surfacetension force remains to be ξ ϕ ∇ ϕ at the gas-liquid interface and smoothly reduces to zero away from it. Asa result, we obtain f s = φξ ϕ ∇ ϕ in Eq.(10). To more clearly illustrate the distribution of the surface tensionforce, we consider a gas bubble at the center of a unit domain with a radius 0 .
2, the bottom half of which isin contact with the solid phase while the upper half is in contact with the liquid phase. Fig.1 schematicallyshows the magnitude of the surface tension forces. We will further investigate the surface tension force inSection 4.1.5. It should be noted that the thermo-capillary effect has not been considered in the presentstudy and the surface tension σ is treated as a constant.To model the no-slip boundary condition at both the gas-solid and liquid-solid interfaces, the velocityinside the solid phase needs to be the given value u S . We add a drag force, formulated as f d = A d ( u S − u ), tothe momentum equation Eq.(10), where the drag coefficient A d can be considered as the Lagrange multiplierthat enforces u = u S inside the solid phase. A d should depend on the volume fraction of the solid phase insuch a way that it is zero away from the solid phase and increases to a large enough value that overwhelms theinertia 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 [89]. Voller and Prakash [90] modified the well-known Carman-Kozeny equation [12], by adding a small constant, denoted by e d here, at the denominator toavoid division by zero, and obtained A d = C d (1 − φ ) φ + e d in Ω M , where C d should be a large number and φ is thevolume fraction of the liquid phase of “ M ” in Ω M as a reminder. Such a definition of A d has been popularlyused, e.g., in [9, 80, 71, 64, 88]. It should be noted that Voller and Prakash [90] did not consider the gasphase and formulated A d in terms of the volume fraction of the liquid phase. Therefore, their formulationis only applicable in Ω M . To obtain A d in Ω, we directly use the volume fraction of the solid phase in Ωas the dependent variable of A d , i.e., A d = C d α S (1 − α S ) + e d , as shown in Eq.(10). Another popular option ofenforcing u = 0 is to assign a large viscosity inside the solid phase, e.g., in [89, 97, 96]. Although it lookssimpler, our numerical tests in Section 4.1.5 show that this is not an effective choice and adding a drag forceis recommended. 8 .1.6 The energy equation The enthalpy of the multiphase mixture is (cid:82) Ω (cid:0) ( ρC p ) T + ρ LM L ( ϕφ ) (cid:1) d Ω, and only the heat conductivity isconsidered in the present study. As a result, we obtain the following energy equation: ∂ (( ρC p ) T ) ∂t + ∇ · ( u ( ρC p ) T ) + ρ LM L (cid:20) ∂ ( ϕφ ) ∂t + ∇ · ( u ϕφ ) (cid:21) = ∇ · ( κ ∇ T ) + Q T . (11)Here, ( ρC p ) and κ denote the volumetric heat and the heat conductivity, respectively, and they are computedfrom Eq.(5). Q T 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 theenthalpy. Therefore, without the phase change, i.e., φ ≡ T > T M or φ ≡ T < T M , these termsshould disappear. To achieve this goal, ∇ · ( u ϕφ ) in Eq.(11) needs to be replaced by ∇ · ( m ϕ φ ) or ∇ · ( m ϕφ ),due to the consistency of volume fraction conservation . As a result, the terms in the bracket in Eq.(11) afterthe 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 T ≡ T in Eq.(11), we obtain T (cid:104) ∂ ( ρC p ) ∂t + ∇ · ( u ( ρC p )) (cid:105) on the left-hand side and 0 on the right-hand side. However, (cid:104) ∂ ( ρC p ) ∂t + ∇ · ( u ( ρC p )) (cid:105) is not zero with ( ρC p )computed from Eq.(5). Therefore, the isothermal solution, i.e., T ≡ T , is not admissible by Eq.(11).Following the derivation from the consistency of mass conservation in Section 2.1.4 and replacing ρ with( ρC p ), we can show that (cid:104) ∂ ( ρC p ) ∂t + ∇ · m ( ρC p ) (cid:105) is zero, where m ( ρC p ) is the flux of volumetric heat definedin 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), ∇ · ( u ( ρC p ) T ) needs to be replaced by ∇ · ( m ( ρC p ) T ). Combining the abovemodifications to Eq.(11), we obtain the consistent energy equation: ∂ (( ρC p ) T ) ∂t + ∇ · ( m ( ρC p ) T ) + ρ LM L (cid:18) ∂ ( ϕφ ) ∂t + ∇ · m ϕφ (cid:19) = ∇ · ( κ ∇ T ) + Q T . (12) Eq.(1), Eq.(4), Eq.(5), Eq.(8), Eq.(10), and Eq.(12) complete the consistent and conservative model forthermo-gas-liquid-solid flows including liquid-solid phase change, and the proposed model honors manyphysical properties. From Eq.(9), Eq.(10), and Eq.(12), it is obvious that the mass and enthalpy of themultiphase mixture are conserved, and the momentum (neglecting the gravity) is conserved without theappearance of the solid phase, by noticing that f s = ξ ϕ ∇ ϕ is equivalent to ∇ · ( − λ ϕ ∇ ϕ ⊗ ∇ ϕ ), see [44, 43,79, 38], and f d = , in this circumstance. When there is the solid phase, the momentum is not necessarilyconserved due to the no-slip boundary condition at the solid boundary. The Galilean invariance is alsosatisfied by the proposed model. The proof is straightforward, using the Galilean transformation, andexamples are available in [41, 37]. The subtle part of the proof is related to the left-hand side of themomentum equation, and we need to emphasize that the consistency conditions are playing a critical rolethere.More importantly, the proposed model is reduction consistent with (i) the isothermal consistent andconservative Phase-Field method for two-phase incompressible flows in [41] when the solid phase is absentand the initial homogeneous temperature is larger than the melting temperature, (ii) the isothermal fictitiousdomain Brinkman penalization (FD/BP) method for fluid-structure interaction in [4, 6] when the liquid phaseis absent and the initial homogeneous temperature is lower than the melting temperature, and (iii) the Phase-Field model of solidification in [7] when the gas phase and flow are absent and the material properties of theliquid 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 conser-vative Phase-Field method for two-phase incompressible flows in [41].Proof.
Given φ = 1 and T = T > T M at t = 0, as already analyzed in Section 2.1.2 and Section 2.1.6,we obtain φ = 1 and therefore m ϕφ = m ϕ and ξ φ = 0 from Eq.(4), and T = T from Eq.(12), at ∀ t > β = β G + ( β LM − β G ) ϕ and m β = β G u +( β LM − β G ) m ϕ , showing that the contribution of β SM disappears. Finally in Eq.(10), f s becomes ξ ϕ ∇ ϕ , and f d = due to α S = ( ϕ − ϕφ ) = 0. Therefore, the temperature remains its homogeneous initial value, andthe simplified system from the proposed model in Section 2.1 with the given condition is equivalent to theconsistent and conservative Phase-Field method for two-phase incompressible flows in [41], by noticing that˜ ϕ = ϕ is the order parameter of the Cahn-Hilliard equation in [41]. Theorem 2.2.
The proposed model in Section 2.1 is consistent with the isothermal fictitious domainBrinkman penalization (FD/BP) method for fluid-structure interaction in [4, 6].Proof.
Given φ = 0 and T = T < T M at t = 0, as already analyzed in Section 2.1.2 and Section 2.1.6,we obtain φ = 0 and therefore m ϕφ = and ξ φ = 0 from Eq.(4), and T = T from Eq.(12), at ∀ t > f s becomes , and f d = A d ( u S − u ) with α S = ϕ .In this case, the gas phase should be understood as an arbitrary incompressible fluid and the solid phaserepresents the fictitious domain, where the material properties are the same as the fluid ones without loss ofgenerality. As a result, we obtain β = β G and m β = β G u from Eq.(5). Therefore, the temperature remainsits initial homogeneous value, and the simplified system from the proposed model in Section 2.1 with thegiven condition is equivalent to the FD/BP method for fluid-structure interaction in [4, 6]. Notice that A d is defined proportional to α S in [4, 6], different from the one in Eq.(10), and the level-set method, insteadof the Cahn-Hilliard equation, is used in [6] for the volume fraction of the fictitious domain (or the solidphase). Theorem 2.3.
The proposed model in Section 2.1 is consistent with the Phase-Field model of solidificationin [7].Proof.
Given ϕ = 1 and u = at t = 0, we have ξ ϕ = 0 and therefore ∂ϕ/∂t = 0 from Eq.(1), which implies ϕ = 1 and m ϕ = u at ∀ t >
0. Further requiring that the material properties of the liquid and solid phasesof “ M ” are identical, we obtain ∇ · u = 0 and ∂ u /∂t = from Eq.(8) and Eq.(10), respectively. Therefore,we obtain u = at ∀ t >
0, and all the convection terms are dropped. Putting all these to the phase changeequation Eq.(4) and the energy equation Eq.(12), they become the same as those in [7], except that p (cid:48) ( φ ) isreplaced with ˜ p (cid:48) ( φ ) in the present study. Remark:
In the proofs of Theorem 2.1, Theorem 2.2, and Theorem 2.3, the given conditions are assumedto be true at t = 0 in the entire domain for convenience. Actually, we only need those conditions to be truelocally 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 wheneverone of the phases is locally absent. The numerical procedure to solve the proposed model in Section 2.1 is introduced in this section. Thedifferential operators are discretized with the conservative finite difference method as those in [41], such thatthe convection terms are discretized by the 5th-order WENO scheme [47], while the divergence, gradient,and Laplacian operators are approximated by the 2nd-order central difference [27]. These discrete operatorshave been carefully validated in various studies [41, 39, 38, 40, 37]. Following the notations in [41], thediscrete operators are denoted by ˜( · ), and time levels are indicated by superscript. We use γ t f n +1 − ˆ f ∆ t todenote the discretization of the time derivative ∂f∂t , and f ∗ ,n +1 to denote the approximation of f n +1 fromprevious time levels. Here, γ t = 1, ˆ f = f n , and f ∗ ,n +1 = f n in the 1st-order case, while they are γ t = 1 . f = 2 f n − . f n − , and f ∗ ,n +1 = 2 f n − f n − in the 2nd order case. Unless otherwise specified, we use the2nd-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 [22, 41], alongwith the consistent and conservative boundedness mapping [39]. Then, the fully-discretized Cahn-Hilliard10quation is recovered and rearranged to be γ t ϕ n +1 − ˆ ϕ ∆ t + ˜ ∇ · ˜ m ϕ = ϕ ∗ ,n +1 ˜ ∇ · u ∗ ,n +1 . (13)Therefore, we obtain ϕ n +1 and the discrete Phase-Field flux ˜ m ϕ after solving the Cahn-Hilliard equation.Second, we proceed to solve the phase change equation Eq.(4). To preserve the consistency of volumefraction conservation on the discrete level, it should be noted that ϕ n +1 and ˜ m ϕ in Eq.(13) are inputs tosolve the phase change equation. The fully-discretized phase change equation is γ t ( ϕ n +1 φ n +1 ) − (cid:91) ( ϕφ )∆ t + ˜ ∇ · ( ˜ m ϕ ˜ φ ∗ ,n +1 ) = − M φ ˜ ξ φ + ϕ ∗ ,n +1 φ ∗ ,n +1 ˜ ∇ · u ∗ ,n +1 , (14)˜ ξ φ = − λ φ ˜ ∇ · ( ϕ n +1 ˜ ∇ φ n +1 ) + λ φ η φ ϕ n +1 ˜ g (cid:48) ( φ n +1 ) + ρ LM LT M ϕ n +1 ˜ p (cid:48) ( φ ∗ ,n +1 )( T M − T ∗ ,n +1 ) , ˜ m ϕφ = ˜ m ϕ ˜ φ ∗ ,n +1 , where ˜ φ represents the WENO reconstruction, ϕ denotes the linear interpolation, and ˜ g (cid:48) ( φ n +1 ) is linearized g (cid:48) ( φ n +1 ) around φ ∗ ,n +1 from Taylor expansion.After solving Eq.(13) and Eq.(14), ( ρC p ), ˜ m ( ρC p ) , and κ are obtained from Eq.(5), noticing that u ∗ ,n +1 is applied to Eq.(5) due to the consistency of reduction , see [38]. Then, the temperature is updated from thefollowing fully-discretized energy equation: γ t ( ρC p ) n +1 T n +1 − (cid:92) ( ρC p ) T ∆ t + ˜ ∇ · ( ˜ m ( ρC p ) ˜ T ∗ ,n +1 ) + ρ LM L (cid:34) γ t ( ϕ n +1 φ n +1 ) − (cid:91) ( ϕφ )∆ t + ˜ ∇ · ˜ m ϕφ (cid:35) (15)= ˜ ∇ · ( κ n +1 ˜ ∇ T n +1 ) + Q n +1 T , where ˜ T represents the WENO reconstruction and κ denotes the linear interpolation. It should be noted thatthe terms in the bracket in Eq.(15) are identical to the left-hand side of the fully-discretized phase changeequation 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 schemeon a collocated grid [41], which has been carefully analyzed and successfully applied to two- and multi-phase problems [41, 38, 39]. Again, ρ , ˜ m ρ , and µ in the momentum equation are directly computed fromEq.(5) and the surface tension force f n +1 s is computed from its definition in Eq.(10) with the balanced-forcemethod [41, 28], while the drag force f d is treated implicitly. As mentioned in Section 2.1.5, f d should bepredominant over either the inertial or viscous effect inside the solid phase. In other words, from Eq.(10), A d | α S =1 = C d /e d should be much larger than ρ/ ∆ t (inertia) or µ/h (viscous). Here h denotes the grid size.To achieve this goal, we set C d as ( ρ n +1 / ∆ t + µ n +1 /h ), and e d is fixed to be 10 − . Therefore, f d will alwaysbe thousand times larger than both the inertial and viscous effects inside the solid phase, regardless of thenumerical setup or material properties. The divergence of the velocity at the new time level, which appearsin the projection scheme, is determined from Eq.(8), i.e.,˜ ∇ · u n +1 = M φ ( ρ LM − ρ SM ) ρ n +1 ξ n +1 φ , (16)where ξ n +1 φ is obtained from its definition in Eq.(4). On the continuous level, as discussed in Section 2.1.4,with ρ and m ρ from Eq.(5) and ∇ · u in Eq.(8), Eq.(9) is implied. However, this is not necessarily true afterdiscretization when the phase change happens and the densities of the liquid and solid phases are not thesame. 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 S m u ∗ ,n +1 is added to the momentumequation, where S m is the residual of the fully-discretized mass conservation equation, i.e., S m = γ t ρ n +1 − ˆ ρ ∆ t + ˜ ∇ · ˜ m ρ . (17)It should be noted that S m only appears on the discrete level due to discretization errors, and it is exactlyzero away from the liquid-solid interface. 11ollowing the above steps, the physical connections among different parts of the proposed model arecorrectly 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 thediscrete level. This will be numerically validated in Section 4.1. In this section, various numerical tests are performed to validate and demonstrate the proposed model inSection 2.1. Then, the predictions from the proposed model are compared to experimental data and othersimulations. 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 beencarefully validated in [41, 39, 38] and therefore not repeated here. Unless otherwise specified, the initialvelocity is zero, and η ϕ = η φ = h and M ϕ λ ϕ = 10 − are set, where h denotes the grid/cell size. We first validate the theorems in Section 2.2 with problems having analytical solutions. Then, the corre-spondence of the mass conservation and the volume change of the phase change material “ M ” is illustrated.Finally, the effectiveness of the proposed surface tension and drag forces is demonstrated. We consider a large-density-ratio advection problem to validate Theorem 2.1 where the solid phase is absentand the temperature is above the melting temperature. The unit domain considered is doubly periodic. Acircular drop, whose density is ρ LM = 10 , is initially at the center of the domain with a radius 0 .
2, surroundedby a gas whose density is ρ G = 1. The specific heats of the phases are 1 × , and the material propertiesof the liquid phase are shared with the solid phase. The viscosity, heat conduction, surface tension, andgravity are neglected. Other parameters are Γ φ = 2 . × − , µ φ = 2 . × − , T M = 2, L = 100, and η ϕ = η φ = 3 h . The solid phase is initially absent, i.e., α S = ϕ − ϕφ = 0, and therefore we have initial φ being 1. The initial velocity and temperature are u = { u , v } = { , } and T = 3 > T M , respectively.The domain is discretized by 128 ×
128 grid cells, and the time step is determined from u ∆ t/h = 0 . α S = ϕ − ϕφ = 0 or φ = 1 at ∀ t >
0, thetemperature remains homogeneous, i.e., T = T at ∀ t >
0, 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 istranslated by the homogeneous velocity. Therefore, there should not be any changes to the shape of the dropand the velocity. From Fig.2 a), we observe that the drop correctly returns to its initial location at t = 1,without any deformation, and the streamlines at t = 1 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) alsoshows that Theorem 2.1 is true due to φ = 1 and T = T at ∀ t > , 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 becomefluctuating, which are observed in [41, 39, 38, 37]. 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 Section2.1.6. A Couette flow problem is solved to validate Theorem 2.2 where the liquid phase is absent and the tem-perature is below the melting temperature. The unit domain considered is periodic along the x axis whileis no-slip along the y axis. Both the top and bottom boundaries are adiabatic but the top one is movingwith a unit horizontal velocity, i.e., u top = 1. The solid phase is at the bottom below y = 0 .
3, while thegas phase fills the rest of the domain. The input parameters are ρ = 1, µ = 0 . C p = 10 , κ = 0,Γ φ = 2 . × − , µ φ = 2 . × − , T M = 2, L = 100, σ = 10 − , and g = . The liquid phase is initially12igure 2: Results of the large-density-ratio advection. a) Streamlines at t = 1 and interface of the drop( ϕ = 0 .
5) at t = 0 and t = 1. Blue arrow lines: streamlines at t = 1, Black solid line: interface at t = 0, Reddashed line: interface at t = 1. The drop returns to its original location at t = 1 without deformation. b) L ∞ norms of ( φ − u − u ), ( v − v ), and ( T − T ) versus time. The velocity preserves its initial value, thesolid phase remains absent, and the temperate equilibrium is maintained exactly, which validate Theorem2.1.absent, i.e., α L = ϕφ = 0, and therefore we have φ = 0. The initial temperature is T = 1 which is lowerthan the melting temperature T M = 2. The domain is discretized by 128 ×
128 grid cells, and the time stepis determined by u top ∆ t/h = 0 . ∂u∂t = ν ∂ u∂y , u = u at y = y , u = u at y = y , u | t =0 = 0 , ν = 0 . , (18)where u = 0, u = u top = 1, y = 0 .
3, and y = 1. The exact solution of Eq.(18) is u E = (cid:80) ∞ α =1 2( u cos( απ ) − u ) απ exp (cid:20) − ν (cid:16) απy − y (cid:17) t (cid:21) sin (cid:104) απy − y ( y − y ) (cid:105) + u − u y − y ( y − y ) + u , y (cid:54) y (cid:54) y , , else , (19)derived from separation of variables. Theorem 2.2 implies that the liquid phase remains absent, i.e., α L = ϕφ = 0 or φ = 0 at ∀ t >
0, the temperature remains homogeneous, i.e., T = T at ∀ t >
0, and the solutionin Eq.(19) is produced (or approximated) by the FD/BP FSI formulation, with the present setup. Expectedresults are obtained and shown in Fig.3. First in Fig.3 a), the profile of the solid fraction ( α S = ϕ − ϕφ = ϕ )at t = 1 overlaps the one at t = 0, representing that the solid phase is stationary. Moreover, the profiles ofthe solution from the proposed model are indistinguishable from the exact solution in Eq.(19), noticing thatthe summation in Eq.(19) is from α = 1 to α = 10 . Theorem 2.2 is true, as shown in Fig.3 b) that φ = 0and T = T at ∀ t >
0. In addition, the unidirectional condition, which is required to obtain Eq.(18), is alsodemonstrated, as shown in Fig.3 b) that v = 0 at ∀ t > The Stefan problem is performed to validate Theorem 2.3 where both the gas phase and the flow are absentand the liquid and solid phases have matched material properties (except the thermal conductivities). Thesetup in [45] is followed. The unit domain is periodic along the x direction. Both the top and bottom13igure 3: Results of the Couette flow. a) Profile of ϕ at t = 0 and t = 1, and profiles of u E and u at t = 0, t = 0 . t = 0 . t = 0 .
6, and t = 1. Here, u E is the exact solution of the Couette flow in Eq.(19), and thenumerical predictions and analytical solutions are overlapped. b) L ∞ norms of φ , v , and ( T − T ) versustime. The liquid phase remains absent, the temperature equilibrium is preserved, and the unidirectional flowcondition is produced, which validate Theorem 2.2.boundaries are free-slip and adiabatic. The material properties are: ρ = 1, µ = 1, C p = 1, κ LM = 0 . κ SM = 1, T M = 1, L = 0 .
53, and g = . Other parameters are η φ = (cid:112) a φ ζ φ , M φ = 1 / ( ν φ ζ φ L ), and λ φ = 1 / ( ν φ M φ ), where a φ = 0 . ν φ = 1, and ζ φ = 0 . y = s = 0 .
2, above which there is the solid phase having a temperature T S = 1 . T L = 1 .
53. The domain is discretized by5 × t = 5 × − .The above setup is to produce the following Stefan problem [45] (based on ∆ T = T − T M ): ∂ ∆ T∂t = ∂∂y (cid:18) κ LM ∂ ∆ T∂y (cid:19) y ∈ ( −∞ , s ( t )) , ∂ ∆ T∂t = ∂∂y (cid:18) κ SM ∂ ∆ T∂y (cid:19) y ∈ ( s ( t ) , + ∞ ) , (20)∆ T = T L − T M at y → −∞ , ∆ T = T S − T M at y → + ∞ , ∆ T = 0 , at y = s ( t ) ,L dsdt = κ SM ∂ ∆ T∂y | y = s + − κ LM ∂ ∆ T∂y | y = s − , where s ( t ) is the location of the liquid-solid interface. Eq.(20) has an analytical self-similar solution [55, 45]: s E ( t ) = s + 2 α √ t, (21) T E − T M = ( T L − T M ) erfc (cid:32) y − s √ κLM t (cid:33) − erfc (cid:32) α √ κLM (cid:33) − erfc (cid:32) α √ κLM (cid:33) y < s ( t ) , ( T S − T M ) − erfc (cid:32) y − s √ κSM t (cid:33) erfc (cid:32) α √ κSM (cid:33) y > s ( t ) ,α = (cid:113) κ SM √ πL T SM − T M erfc( α/ (cid:113) κ SM ) exp (cid:18) − α κ SM (cid:19) + (cid:113) κ LM √ πL T LM − T M − erfc (cid:18) α/ (cid:113) κ LM (cid:19) exp (cid:18) − α κ LM (cid:19) . Numerical results are compared to the exact solution Eq.(21), and shown in Fig.4. The interface location is14igure 4: Results of the Stefan problem. a) Profiles of the temperature at different time. Here, T E is theexact solution of the Stefan problem in Eq.(21). b) Location of the liquid-solid interface versus time. Thenumerical prediction agrees well with the analytical solution, which validates Theorem 2.3.specified as the 0 . ϕφ ) from the numerical results. Both the temperature and interface locationagree with the exact solution very well. Minor discrepancy is observed near the domain boundary at t = 0 . L ∞ norms of ( ϕ − u , and v are on the orders of 10 − , 10 − , and 10 − , respectively, at the end of the simulation, which demonstratesTheorem 2.3. Here, we consider the effect of the mass conservation, which leads to the non-divergence-free velocity, i.e.,Eq.(8), when the phase change happens and the liquid and solid phases have different densities. Such aneffect has been overlooked by many existing models. The unit domain considered is periodic along the x axis. The bottom boundary is no-slip and has a fixed temperature T bottom = 5, while the top one isan outflow boundary having a fixed pressure P top = 0 and zero heat flux. At the bottom of the domainbelow y = 0 . ρ LM = 1000, µ LM = 1 × − , ( C p ) LM = 1,and κ LM = 200. Floating on the liquid phase is the solid phase, whose material properties are ρ SM = 900, µ SM = 1 × − , ( C p ) SM = 1, and κ SM = 300. Above y = 0 . ρ G = 1, µ G = 2 × − , ( C p ) G = 0 .
1, and κ G = 100. Other parameters are T M = 1, L = 20, σ = 0 . ρ LM L Γ φ /T M = 10, M φ = 5 × − , and g = (0 , − . T = 0 .
5. The domain isdiscretized by 128 ×
128 grid cells, and the time step is ∆ t = 10 − .The materials are heated by the bottom wall whose temperature is higher than the melting temperature.As a result, the solid phase will melt and finally disappear. Since the liquid density is 10% larger than thesolid phase, the final volume of the phase change material should be smaller than its initial value, in orderto honor the mass conservation. Results are shown in Fig.5 and match the expectation. From Fig.5 a), theliquid-solid interface is moving upward while at the same time the gas-solid interface is moving downward.At the beginning, the volume (area) of “ M ”, including its liquid and solid phases, is 0 .
6. At the end of thesimulation, there is only the liquid phase of “ M ” in the domain, and the gas-liquid interface stays horizontallybelow y = 0 .
6, indicating that the volume of “ M ” is smaller than its initial value.Quantitative data are reported in Fig.5 b) where the displacements of the liquid-solid and gas-solidinterfaces versus time are plotted. The gas-solid interface is defined as the 0 . ϕ , while theliquid-solid interface is the 0 . ϕφ ). We observe that the liquid-solid interface actually movesdownward at the very beginning because the initial temperature is below the melting temperature. Asa result, solidification happens in that period. As the materials are heated from the bottom wall, thesolid melts, leading to the rise of the liquid-solid interface but fall of the gas-solid interface, as expected.The melting process ends before t = 4. We estimate the mass of the phase change material “ M ” simply15y [ ρ LM s LS + ρ SM ( s GS − s LS )], where s LS and s GS denote the locations of the liquid-solid and gas-solidinterfaces, respectively, illustrated in the first snapshot in Fig.5 a). As plotted in Fig.5 b), the change of[ s GS + ( ρ LM /ρ SM − s LS ] is negligible, which implies that the movements of the interfaces are constrainedby the mass conservation. One can expect a more obvious displacement of the gas-solid interface, inducedby the phase change, to appear if the density difference of the liquid and solid phases of the phase changematerial is larger than the present setup.Although the proposed model strictly satisfies the mass conservation, i.e., Eq.(9), the present schemedoes not always do, as discussed in Section 3. Fig.5 c) shows the relative changes of the total mass ( (cid:82) Ω ρd Ω)and the mass of “ M ” ( (cid:82) Ω ρ M d Ω) versus time. Here, ρ M = ρ LM ( ϕφ ) + ρ SM ( ϕ − ϕφ ) and the integral iscomputed from the mid-point rule. It should be noted that the total mass in this case is not conserved, andits change is related to the volume change of “ M ” during the phase change. The initial decrease of the totalmass corresponds to the solidification process, where the volume of “ M ” expands and therefore the gas issqueezed out. As melting occurs, the total mass increases because the volume of “ M ” reduces, and the gasmoves into the domain. When melting is completed, the total mass stops changing as well. On the otherhand, the mass of “ M ” should be conserved even though the phase change happens, while it is not exactlytrue due to numerical errors. Nonetheless, its relative change is very small, on the order of 10 − , whichis satisfactory. It should be noted that as long as the velocity is divergence-free, i.e., in the present workthe phase change is absent or the densities of the liquid and solid phases are the same, the present schemeexactly conserves the mass of “ M ” as well as “ G ” on the discrete level, see [41, 39].In summary, the proposed model automatically and physically captures the volume change induced bythe phase change, and therefore the mass conservation, thanks to the consistency of mass conservation . Thisphysical behavior is not correctly captured in many existing models, where the velocity is assumed to bedivergence-free. Here, we demonstrate the performances of the surface tension force f s , which models the surface tension atthe gas-liquid interface, and the drag force f d , which enforces zero velocity in the solid phase and thereforethe no-slip condition at the solid boundary, in the momentum equation Eq.(10). Unless otherwise specifiedin this section, the following setup is employed. The material properties are ρ LM = 2 . × kg / m , µ LM = 1 . × − Pa · s, ( C p ) LM = 1 . × J / (K · kg), κ LM = 91W / (m · K), ρ SM = 2 . × kg / m , µ SM =1 . × − Pa · s, ( C p ) SM = 0 . × J / (K · kg), κ SM = 211W / (m · K), ρ G = 0 . / m , µ G = 4 × − Pa · s,( C p ) G = 1 . × J / (K · kg), κ G = 61 × − W / (m · K), Γ φ = 1 . × − m · K, µ φ = 1 . × − m / (s · K), T M = 933 . L = 3 . × J / kg, σ = 0 . / m, and g = (0 , − . / m , a length scale 0 . / s , and atemperature scale 933 . × T bottom = 0 .
5. The center of a circular bubble havinga radius 0 .
125 is at (0 . , . y = 0 .
1, while the liquid phase fills the rest ofthe domain. The initial temperature is 0 . . ×
128 grid cells, and the timestep is ∆ t = 1 × − .In the first two cases, the drag force is zero, i.e., f d = , while the solid viscosity becomes 1Pa · s, whichis about 1000 times larger than the liquid phase. In case 1, we employ the proposed surface tension force f s = φξ ϕ ∇ ϕ in Eq.(10), while it is f s = ξ ϕ ∇ ϕ in case 2. Results are shown in Fig.6, and the differencebetween case 1 ( f s = φξ ϕ ∇ ϕ ) and case 2 ( f s = ξ ϕ ∇ ϕ ) is obvious. In case 1, the surface tension force onlyacts at the upper part of the bubble, contacting the liquid phase, which is desirable. As a result, the bottompart of the bubble is easier to be deformed, while the upper part tends to be flattened. On the other handin case 2, the surface tension force acts on the entire bubble interface, no matter whether the bubble iscontacting the liquid or solid phase. Consequently, the bubble remains circular even after it is contactedby the solid phase at its bottom part. The above analysis is demonstrated in Fig.7, where the magnitudeof the surface forces, i.e., | f s | , at t = 0 .
20 in cases 1 and 2 is shown. It can be learned from Fig.6 thatthe surface tension force can be influential to the results, and Fig.7 demonstrates that the proposed surfacetension force, i.e., f s = φξ ϕ ∇ ϕ in case 1, is the one that should be chosen. An alternative formulation, i.e., f s = ξ ϕ ∇ ( ϕφ ), has also been tested, and it produced an unstable solution. Another issue, observed in Fig.6,16igure 5: Results of the mass conservation and volume change. a) Snapshots of the phases at t = 0 . t = 2 .
5, and t = 5 .
0, and schematic of the locations of the liquid-solid ( s LS ) and gas-solid ( s GS ) interfaces.White: the gas phase, Orange: the liquid phase, Blue: the solid phase. b) Displacements of the gas-solidand liquid-solid interfaces versus time. The displacements of the interfaces are constrained by the massconservation, quantified by ∆[ s GS + ( ρ LM /ρ SM − s LS ] = 0. c) Relative changes of the total mass and themass of the phase change material “ M ” versus time. The total mass is changing because the gas moves inand out of the domain following the volume change of the phase change material. The mass of the phasechange material changes slightly (less than 0 . . × times more than the gas phase. In the ideal situation, the viscous force is infinite insidethe solid phase, which in turn enforces zero velocity gradient there. In other words, after discretization,˜ µ SM /h >> ˜ ρ SM / ∆ t should be true, where ˜ µ and ˜ ρ are non-dimensionalized viscosity and density, and inthis specific case, µ SM >> . · s. We again tried µ SM = 1000Pa · s and it quickly became unstable.Therefore, increasing the solid viscosity is not an effective way to enforce zero velocity in the solid phase.This is the reason the drag force f d is introduced in the proposed model.Next, the drag force f d is activated. In case 3, we apply the formulation of A d in Eq.(10), while in case 4,an alternative definition of A d , i.e., A d = ϕC d (1 − φ ) φ + e d , is considered. The alternative A d can be easily derivedfrom the drag force model in Ω M proposed by Voller and Prakash [90] and the diffuse domain approach [52].Results are shown in Fig.8. It is obvious that the solid movement is suppressed after comparing Fig.8 toFig.6, and this has also been quantitatively demonstrated in Section 4.1.2. As long as the bubble is “caught”by the solid from the bottom, it stops rising, unlike the one in Fig.6. Although both cases 3 and 4 producesimilar results, one can observe that the bubble in case 4 is less deformed than the one in case 3 using theproposed formulation. This implies that the alternative A d in case 4 has a larger effective region to enforcethe velocity to be zero, while its influence on the overall dynamics is negligible. We conclude that bothchoices of A d are valid, but we keep using the one in Eq.(10) in the present study because it is equivalent to17igure 6: Results of cases 1 and 2 at t = 0 . t = 0 .
05, and t = 0 . t = 0 . t = 0 .
50. Orange: the liquidphase, Blue: the solid phase. Top: case 1 using the proposed surface tension force f s = φξ ϕ ∇ ϕ . Bottom:case 2 using the original surface tension force f s = ξ ϕ ∇ ϕ . The original surface tension force (case 2) alsoappears at the gas-solid interface and therefore the bubble is less deformable than the one in case 1. Thesolid behaves like a fluid although its viscosity is about 1000 times larger than the liquid phase.Figure 7: Magnitude of the surface tension forces | f s | at t = 0 .
20. a) Case 1: f s = φξ ϕ ∇ ϕ , b) Case 2: f s = ξ ϕ ∇ ϕ . The proposed surface tension force in case 1 appears only at the gas-liquid interface, while theoriginal one in case 2 mistakenly appears at the gas-solid interface.the Carman-Kozeny equation [12]. Here, we compare results from the proposed model to the experimental [29] and numerical [9, 50] results.Details of the setup have been given in [29, 9], and we follow those in the present study. A rectangular cavity,whose width is 8 . . . . .
45K at the right wall, while both the top and bottomwalls are adiabatic. The material properties of gallium are: density 6093kg / m , viscosity 1 . × − Pa · s,specific heat 381 . / (K · kg), thermal conductivity 32W / (m · K), melting temperature 302 . / kg. The gravity is g = (0 , − . / m , and the buoyancy force is computed from the Boussinesqapproximation, i.e., F g = − ρ β ( T − T ) g , where ρ = 6095kg / m , T = 302 . β = 1 . × − / K isthe thermal expansion coefficient. Γ φ and µ φ are chosen to be 1 . × − m · K and 1 . × − m / (s · K). The18igure 8: Results of cases 3 and 4 at t = 0 . t = 0 .
05, and t = 0 . t = 1 . t = 3 .
40. Orange: the liquidphase, Blue: the solid phase. Top: case 3 using A d = C d α S (1 − α S ) + e d . Bottom: case 4 using A d = ϕC d (1 − φ ) φ + e d .Both choices of A d successfully suppress the solid motion and produce similar results, but the one in case 4has a larger effective region resulting in a less-deformed bubble.domain is discretized by 126 ×
96 grid cells. The initial time step is ∆ t = 0 .
05 and adaptively changes to be0 . h/ max( u, v ). The governing equations are non-dimensionalized by a density scale 6095kg / m , a lengthscale 0 . . / s , and a temperature scale 302 . p ( φ ) inEq.(2) is applied, because its derivative is zero in the solid-state. This demonstrates the significance of usingthe proposed ˜ p (cid:48) ( φ ) in Eq.(4) in realistic problems. The more detailed analysis has been provided in Section2.1.2. It is preferable to understand µ φ and Γ φ as tunable parameters of the proposed model, instead of theirphysical meaning, since the practical interface thickness is much larger than the physical value. In practice,we tune µ φ and Γ φ so that the numerical result matches the experimental one at t = 2min, and obtain therest of the results with those parameters. When tuning µ φ and Γ φ , we discover that µ φ controls the speedof the phase change, while Γ φ affects the interface thickness. A larger µ φ gives a larger M φ , and as a resultaccelerates the phase change. We observe an over-compressed interface when Γ φ is too small, while a toolarge Γ φ casts difficulty to initialize the interface. This can be explained by the energy mechanism in thePhase-Field model of solidification Eq.(2). λ φ controls the net effect of the thermodynamical compressionand diffusion that preserve the interface thickness. A too small λ φ , resulting from a small Γ φ , basicallyremoves those compression and diffusion effects and leads to a sharp interface. On the other hand, a large λ φ strengthens those effects, and a larger overheat is therefore needed to drive the order parameter, jumpingfrom one equilibrium state to another across the double-well potential. We suggest µ φ and Γ φ sharing thesame value. We also test the effect of e d in the drag force f d , and little difference is observed when reducing e d from the default value 10 − to 10 − . Here, we consider a complicated case including interactions among the gas, liquid, and solid phases. Thematerial properties of the liquid phase are ρ LM = 2 . × kg / m , µ LM = 1 . × − Pa · s, ( C p ) LM =1 . × J / (K · kg), and κ LM = 91W / (m · K). They are ρ SM = 2 . × kg / m , µ SM = 1 . × − Pa · s,( C p ) SM = 0 . × J / (K · kg), and κ SM = 211W / (m · K) for the solid phase, and ρ G = 0 . × kg / m ,19igure 9: Results of melting of gallium in a rectangular cavity and the liquid-solid interface is presented.Both Brent et al. (1988) [9] and Kim et al. (2011) [50] use the enthalpy-porosity technique but Brent et al.(1988) [9] neglect the convection of the liquid fraction in the energy equation. The present results behavesimilarly to those in Kim et al. (2011) [50] and a stronger melt convection is predicted. Nevertheless, all thenumerical predictions are consistent with the experimental data. µ G = 4 × − Pa · s, ( C p ) G = 1 . × × J / (K · kg), and κ G = 61 × − W / (m · K) for the gas phase.The melting temperature is T M = 933 . L = 3 . × J / kg, the surface tension is σ = 0 . / m, the gravity is g = (0 , − . / s , and the Gibbs-Thomson and linear kinetic coefficients arechosen to be Γ φ = 1 . × − m · K and µ φ = 1 . × − m / (s · K), respectively. The governing equationsare non-dimensionalized by a density scale 1kg / m , a length scale 0 . / s , anda temperature scale 933 .
6K the same as the melting temperature.A unit domain is considered. Both the left and right boundaries are no-slip and adiabatic walls. Thebottom boundary is no-slip with a fixed temperature T bottom = 0 .
5. The top boundary has a fixed pressure P top = 0 and a zero heat flux. The domain is discretized by 128 ×
128 cells, and the time step is ∆ t = 10 − .The initial condition of the phases is illustrated in the first snapshot of Fig.10. Above y = 0 .
75 is the gasphase, while the solid phase is at the bottom below y = 0 .
1. In the middle of the domain is the liquidphase inside which there are three circular gas bubbles. The radii of the bubbles from left to right are0 . . .
1, and their centers are at (0 . , . . , . . , . T bottom , inside the solid phase while it is 1.1 elsewhere. Note that thenon-dimensionalized melting temperature is 1.Results are shown in Fig.10. The gas, liquid, and solid phases are filled by the white, orange, and bluecolors, and the solid phase becomes green when the phase change is finished. The bubbles are moving upward20ue to the buoyancy effect, and, at the same time, the liquid is solidified as its temperature is cooled downby the bottom wall. The motion of the bubbles drives the liquid and produces melt convection. As a result,the liquid below the bubbles is solidified faster than its neighbor, and the liquid-solid interface first “catches”the left bubble then the right one. As the largest bubble at the middle rises, the gas-liquid interface abovestarts to be perturbed, which, in turn, deviates the bubble rising from its vertical line. When the middlebubble merges the gas-liquid interface, a strong capillary wave is produced due to the surface tension. Asthe capillary wave travels back and forth, the liquid-solid interface keeps moving upward. Since the heatconductivity of the gas is much smaller than the liquid or solid, the solidification is slower right above thetwo trapped gas bubbles, and the liquid-solid interface forms a “V” shape there. As the liquid-solid interfacegets closer to the gas-liquid one, the capillary wave is quickly attenuated by the viscosity, due to the zerovelocity of the solid. At the end of the simulation, the liquid is completely solidified with two hollows formedby the right and left bubbles. Here, we consider melting a solid rectangle and solidify it again. The material properties and setup areidentical to those in Section 4.3, except that the thermal conductivity of the gas is κ G = 100W / (m · K)and that the bottom wall becomes adiabatic, and the temperature at the other boundaries is 2 before t = 2then 0 .
5. The initial condition of the phases is illustrated in the first snapshot of Fig.11. A rectangularsolid with width 0 . . .
075 and 0 .
1, and centers are located at (0 . , .
1) and (0 . , . . In the present work, we consider the thermo-gas-liquid-solid flows, where the liquid and solid phases areexperiencing solidification/melting. A novel consistent and conservative Phase-Field model is developed forsuch a kind of problem. The ingredients of the proposed model are the consistent and conservative Phase-Field method for incompressible two-phase flows [41], the fictitious domain Brinkman penalization (FD/BP)method for fluid-structure interactions [4, 6], and the Phase-Field model of solidification in [7]. Thesesuccessful models are physically coupled using the consistency of reduction , consistency of volume fractionconservation , consistency of mass conservation , and consistency of mass and momentum transport . Theseconsistency conditions, which have been successfully applied in isothermal, multiphase, multicomponent,immiscible, and incompressible flows [41, 39, 38, 40, 37], are used in problems having variable temperatureand phase changes for the first time, and are demonstrated to play an essential role in the present study.The Cahn-Hilliard equation Eq.(1) is applied to locate the gas and phase change material. The phase changeequation Eq.(4) is derived from the solidification model in [7] using the diffuse domain approach [52]. Thenthe consistency of volume fraction conservation is applied to admit the fully liquid/solid-state of the phasechange material. The interpolation function in [7] is also modified so that the equilibrium states of the orderparameter in the model depend on the temperature in a physical sense, which resolves the issue of initiatingthe phase change when the phase change material is fully liquid/solid at the beginning. After applyingthe consistency of mass conservation , we not only obtain the consistent mass flux, which appears in themomentum equation following the consistency of mass and momentum transport , but also the divergence ofthe velocity Eq.(8), which quantifies the volume change induced by the solidification/melting. Isothermal (ortemperature equilibrium) solutions are admissible by the proposed energy equation Eq.(12) when the phase21igure 10: Results of rising bubbles with solidification. White: the gas phase, Orange: the liquid phase,Blue: the solid phase, Green: the solid phase when the phase change is finished. From left to right and topto bottom, t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 1 . t = 1 . t = 1 . t = 1 . t = 1 . t = 1 . t = 2 . t = 2 . t = 2 . t = 3 . t = 3 . t = 4 .
00, and t = 5 . consistency of mass conservation and the consistency of volumefraction conservation . To confine the surface tension effect on the gas-liquid interface only, we propose a newcontinuous surface tension force based on the one in [41]. The Carman-Kozeny equation [12] is modified toenforce zero velocity in the solid phase. These two additional forces are added to the momentum equationEq.(10). The proposed model defines the volume fractions of the gas, liquid, and solid phases unambiguously,and the volume change due to solidification/melting is included. The mass and energy conservation is alwaystrue, while the momentum conservation is honored if the solid phase is absent due to the no-slip condition atthe solid boundary. The proposed model also satisfies the Galilean invariance. Moreover, we show Theorem22igure 11: Results of rising bubbles with solidification. White: the gas phase, Orange: the liquid phase,Blue: the solid phase, Green: the solid phase when the phase change is finished. From left to right and topto bottom, t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 0 . t = 1 . t = 1 . t = 1 . t = 1 . t = 1 . t = 1 . t = 2 . t = 2 . t = 2 . t = 2 . t = 2 . t = 2 . t = 2 . t = 3 . t = 3 . t = 3 . t = 3 .
75, and t = 4 . C d in the drag force model.After validating the proposed model, a comparison to experimental and other numerical data is conducted,and the proposed model produces results that agree well with those data. We discover that the Gibbs-Thomson and linear kinetic coefficients, i.e., Γ φ and µ φ , need to be carefully selected to obtain a quantitativeagreement. µ φ and Γ φ is positively correlated to the speed of phase change and the interface thickness,respectively, which can be explained by the energy mechanism of the Phase-Field model of solidificationin [7]. In practice, calibration may be needed to determine µ φ , and we suggest that µ φ and Γ φ take thesame value. Finally, two challenging problems, including a wide range of material properties and stronginteractions among different phases, are set up and successfully solved, which illustrates the capability of themodel.The present study proposes a practical framework to incorporate the solidification/melting of a purematerial into liquid-gas flows. This method can be extended to include more complicated physics, e.g.,the thermo-capillary effect, anisotropy or dendritic growth, and solute transport during the solidification.Adaptive grid refinement and time stepping are attractive directions to improve numerical simulations, whilelimiting computational cost. However, the physical properties of the model need to be preserved on thediscrete level to avoid unphysical behaviors. This is a non-trivial problem and deserves further investigation. Acknowledgments
A.M. Ardekani would like to acknowledge the financial support from the National Science Foundation (CBET-1705371). This work used the Extreme Science and Engineering Discovery Environment (XSEDE) [84],which is supported by the National Science Foundation grant number ACI-1548562 through allocation TG-CTS180066 and TG-CTS190041. G. Lin would like to acknowledge the support from National Science Foun-dation (DMS-1555072 and DMS-1736364, CMMI-1634832 and CMMI-1560834), and Brookhaven NationalLaboratory Subcontract 382247, ARO/MURI grant W911NF-15-1-0562, and U.S. Department of Energy(DOE) Office of Science Advanced Scientific Computing Research program DE-SC0021142.24 eferences [1] H. Abels, H. Garcke, and G. Grun. Thermodynamically consistent, frame indifferent diffuse interfacemodels for incompressible two-phase flows with different densities.
Mathematical Models and Methodsin Applied Sciences , 22:1150013, 2012.[2] S.M. Allen and J.W. Cahn. A microscopic theory for antiphase boundary motion and its applicationto antiphase domain coarsening.
Acta Metallurgica , 27:1085–1095, 1979.[3] D.M. Anderson, G.B. McFadden, and A.A. Wheeler. Diffuse-interface methods in fluid mechanics.
Annu. Rev. Fluid Mech. , 30:139–165, 1998.[4] P. Angot, C.-H. Bruneau, and P. Fabrie. A penalization method to take into account obstacles inincompressible viscous flows.
Numerische Mathematik , 81(4):497–520, 1999.[5] C. Beckermann, H.-J. Diepers, I. Steinbach, A. Karma, and X. Tong. Modeling melt convection inphase-field simulations of solidification.
Journal of Computational Physics , 154(2):468–496, 1999.[6] M. Bergmann and A. Iollo. Modeling and simulation of fish-like swimming.
Journal of ComputationalPhysics , 230(2):329–348, 2011.[7] W.J. Boettinger, J.A. Warren, C. Beckermann, and A. Karma. Phase-field simulation of solidification.
Annual review of materials research , 32(1):163–194, 2002.[8] J.U. Brackbill, D.B. Kothe, and C. Zemach. A continuum method for modeling surface tension.
J.Comput. Phys. , 100:335–354, 1992.[9] A.D. Brent, V.R. Voller, and K.T.J. Reid. Enthalpy-porosity technique for modeling convection-diffusion phase change: application to the melting of a pure metal.
Numerical Heat Transfer, Part AApplications , 13(3):297–318, 1988.[10] M. Bussmann, D.B. Kothe, and J.M. Sicilian. Modeling high density ratio incompressible interfacialflows. In
Proceedings of ADME Fluid Engineering Division Summer Meeting , page 31125, 2002.[11] J.W. Cahn and J.E. Hilliard. Free energy of a nonuniform system, i interfacial free energy.
J. Chem.Phys. , 28:258–267, 1958.[12] P.C. Carman. Fluid flow through granular beds.
Chemical Engineering Research and Design , 75:S32–S48, 1997.[13] C. Chan, J. Mazumder, and M.M. Chen. A two-dimensional transient model for convection in lasermelted pool.
Metallurgical Transactions A , 15(12):2175–2184, 1984.[14] C. Chen and X. Yang. Efficient numerical scheme for a dendritic solidification phase field model withmelt convection.
Journal of Computational Physics , 388:41–62, 2019.[15] L.-Q. Chen. Phase-field models for microstructure evolution.
Annual review of materials research ,32(1):113–140, 2002.[16] V.L. Chenadec and H. Pitsch. A monotonicity preserving conservative sharp interface flow solver forhigh density ratio two-phase flows.
J. Comput. Phys. , 249:185–203, 2013.[17] R. Chiodi and O. Desjardins. A reformulation of the conservative level set reinitialization equation foraccurate and robust simulation of complex multiphase flows.
J. Comput. Phys. , 343:186–200, 2017.[18] P-H Chiu and Y-T Lin. A conservative phase-field method for solving incompressible two-phase flows.
J. Comput. Phys. , 230:185–204, 2011.[19] J.A. Dantzig and M. Rappaz.
Solidification: -Revised & Expanded . EPFL press, 2016.2520] N.S. Dhaidan and J.M. Khodadadi. Melting and convection of phase change materials in differentshape containers: A review.
Renewable and Sustainable Energy Reviews , 43:449–477, 2015.[21] H. Ding, P.D.M. Spelt, and C. Shu. Diffuse interface model for incompressible two-phase flows withlarge density ratios.
J. Comput. Phys. , 226:2078–2095, 2007.[22] S. Dong and J. Shen. A time-stepping scheme involving constant coefficient matrices for phase-field sim-ulations of two-phase incompressible flows with large density ratios.
J. Comput. Phys. , 231:5788–5804,2012.[23] Y. Dutil, D.R. Rousse, N.B. Salah, S. Lassue, and L. Zalewski. A review on phase-change materials:Mathematical modeling and simulations.
Renewable and sustainable Energy reviews , 15(1):112–130,2011.[24] B. Echebarria, R. Folch, A. Karma, and M. Plapp. Quantitative phase-field model of alloy solidification.
Physical review E , 70(6):061604, 2004.[25] R.P. Fedkiw, T. Aslam, B. Merriman, and S. Osher. A non-oscillatory eulerian approach to interfacesin multimaterial flows (the ghost fluid method).
J. Comput. Phys. , 152:457–492, 1999.[26] J.J. Feng, C. Liu, J. Shen, and P. Yue. An energetic variational formulation with phase field methodsfor interfacial dynamics of complex fluids: advantages and challenges. In
Modeling of soft matter , pages1–26. Springer, 2005.[27] J.H. Ferziger and M. Peric.
Computational Methods for Fluid Dynamics . Springer Berlin / Heidelberg,third rev. edition, 2001.[28] M.M. Francois, J.S. Cummins, E.D. Dendy, D.B. Kothe, M.J. Sicilian, and W.W. Williams. A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume trackingframework.
J. Comput. Phys. , 213:141–173, 2006.[29] C. Gau and R. Viskanta. Melting and solidification of a pure metal on a vertical wall.
Journal of HeatTransfer , 108:174–181, 1986.[30] F. Gibou, R. Fedkiw, and S. Osher. A review of level-set methods and some recent applications.
J.Comput. Phys. , 353:82–109, 2018.[31] D. Gueyffier, J. Li, A. Nadim, R. Scardovelli, and S. Zaleski. Volume-of-fluid interface tracking withsmoothed surface stress methods for three-dimensional flows.
J. Comput. Phys. , 152:423–456, 1999.[32] Q. He, H. Xia, J. Liu, X. Ao, and S. Lin. Modeling and numerical studies of selective laser melting:Multiphase flow, solidification and heat transfer.
Materials & Design , 196:109115, 2020.[33] C.W. Hirt and B.D. Nichols. Volume of fluid (vof) method for the dynamics of free boundaries.
J.Comput. Phys. , 39:201–225, 1981.[34] S.F. Hosseinizadeh, F.L. Tan, and S.M. Moosania. Experimental and numerical studies on performanceof pcm-based heat sink with different configurations of internal fins.
Applied Thermal Engineering ,31(17-18):3827–3838, 2011.[35] H. Hu and S.A. Argyropoulos. Mathematical modelling of solidification and melting: a review.
Mod-elling and Simulation in Materials Science and Engineering , 4(4):371, 1996.[36] T.-H. Huang, T.-H. Huang, Y.-S. Lin, C.-H. Chang, P.-Y. Chen, S.-W. Chang, and C.-S. Chen.Phase-field modeling of microstructural evolution by freeze-casting.
Advanced Engineering Materi-als , 20(3):1700343, 2018.[37] Z. Huang, G. Lin, and A.M. Ardekani. A consistent and conservative model and its scheme for n-phase-m-component incompressible flows. arXiv:2101.04252 , 2020.2638] Z. Huang, G. Lin, and A.M. Ardekani. A consistent and conservative phase-field method for multiphaseincompressible flows. arXiv:2010.01099 , 2020.[39] Z. Huang, G. Lin, and A.M. Ardekani. Consistent and conservative scheme for incompressible two-phase flows using the conservative allen-cahn model.
J. Comput. Phys. , 420:109718, 2020.[40] Z. Huang, G. Lin, and A.M. Ardekani. A consistent and conservative volume distribution algorithmand its applications to multiphase flows using phase-field models. arXiv:2010.01738 , 2020.[41] Z. Huang, G. Lin, and A.M. Ardekani. Consistent, essentially conservative and balanced-force phase-field method to model incompressible two-phase flows.
J. Comput. Phys. , 406:109192, 2020.[42] Satoshi Ii, Kazuyasu Sugiyama, Shintaro Takeuchi, Shu Takagi, Yoichiro Matsumoto, and Feng Xiao.An interface capturing method with a continuous function: The thinc method with multi-dimensionalreconstruction.
J. Comput. Phys. , 231(5):2328–2358, 2012.[43] D. Jacqmin. Calculation of two-phase navier-stokes flows using phase-field modeling.
J. Comput. Phys. ,155:96–127, 1999.[44] D. Jamet, D. Torres, and J.U. Brackbill. On the theory and computation of surface tension: The elim-ination of parasitic currents through energy conservation in the second-gradient method.
J. Comput.Phys. , 182:262–276, 2002.[45] E. Javierre, C. Vuik, F.J. Vermolen, and S. Van der Zwaag. A comparison of numerical models forone-dimensional stefan problems.
Journal of Computational and Applied Mathematics , 192(2):445–459,2006.[46] Y. Ji, L. Chen, and L.-Q. Chen. Understanding microstructure evolution during additive manufacturingof metallic alloys using phase-field modeling. In
Thermo-Mechanical Modeling of Additive Manufac-turing , pages 93–116. Elsevier, 2018.[47] G-S Jiang and C-W Shu. Efficient implementation of weighted eno schemes.
J. Comput. Phys. ,126:202–228, 1996.[48] S.G. Kim and W.T. Kim. Phase-field modeling of solidification. In
Handbook of materials modeling ,pages 2105–2116. Springer, 2005.[49] Y. Kim, A. Hossain, and Y. Nakamura. Numerical study of melting of a phase change material (pcm)enhanced by deformation of a liquid–gas interface.
International Journal of Heat and Mass Transfer ,63:101–112, 2013.[50] Y.K. Kim, A. Hossain, S. Kim, and Y. Nakamura. A numerical study on time-dependent melting anddeformation processes of phase change material (pcm) induced by localized thermal input.
Two phaseflow phase change and numerical modeling , 23, 2011.[51] B. Lalanne, L.R. Villegas, S. Tanguy, and F. Risso. On the computation of viscous terms for incom-pressible two-phase flows with level set/ghost fluid method.
J. Comput. Phys. , 301:289–307, 2015.[52] X. Li, J. Lowengrub, A. Ratz, and A. Voigt. Solving pdes in complex geometries: A diffuse domainapproach.
Commun. Math. Sci. , 1:81–107, 2009.[53] S. Lin, Z. Gan, J. Yan, and G.J. Wagner. A conservative level set method on unstructured meshesfor modeling multiphase thermo-fluid flow in additive manufacturing processes.
Computer Methods inApplied Mechanics and Engineering , 372:113348, 2020.[54] L.-X. Lu, N. Sridhar, and Y.-W. Zhang. Phase field simulation of powder bed-based additive manu-facturing.
Acta Materialia , 144:801–809, 2018.[55] James M.H.
One-dimensional Stefan problems: an introduction , volume 31. Longman Sc & Tech, 1987.2756] S. Mirjalili, S. Jain, and Dodd M.S. Interface-capturing methods for two-phase flos: An overview andrecent developments.
Center for Turbulence Research Annual Research Briefs , pages 117–135, 2017.[57] Abu-Al-Saud M.O., S. Popinet, and H.A. Tchelepi. A conservative and well-balanced surface tensionmodel.
J. Comput. Phys. , 371:896–931, 2018.[58] N. Nangia, E.G. Boyce, N.A. Patankar, and A.P.S. Bhalla. A robust incompressible navier-stokessolver for high density ratio multiphase flows.
J. Comput. Phys. , 390:548–594, 2019.[59] B. Nestler, A.A. Wheeler, L. Ratke, and C. St¨ocker. Phase-field model for solidification of a monotecticalloy with convection.
Physica D: Nonlinear Phenomena , 141(1-2):133–154, 2000.[60] E. Olsson and G. Kreiss. A conservative level set method for two phase flow.
J. Comput. Phys. ,210:225–246, 2005.[61] E. Olsson, G. Kreiss, and S. Zahedi. A conservative level set method for two phase flow ii.
J. Comput.Phys. , 225:785–807, 2007.[62] S. Osher and A.J. Sethian. Fronts propagating with curvature-dependent speed: Algorithms based onhamilton-jacobi formulations.
J. Comput. Phys. , 79:12–49, 1988.[63] M. Owkes and O. Desjardins. A mass and momentum conserving unsplit semi-lagrangian frameworkfor simulating multiphase flows.
J. Comput. Phys. , 332:21–46, 2017.[64] C. Panwisawas, C. Qiu, M.J. Anderson, Y. Sovani, R.P. Turner, M.M. Attallah, J.W. Brooks, and H.C.Basoalto. Mesoscale modelling of selective laser melting: Thermal fluid dynamics and microstructuralevolution.
Computational Materials Science , 126:479–490, 2017.[65] W. Pitscheneder, T. DebRoy, K. Mundra, and R. Ebner. Role of sulfur and processing variables on thetemporal evolution of weld pool geometry during multikilowatt laser beam welding of steels.
WeldingJournal (Miami, Fla) , 75(3):71s–80s, 1996.[66] S. Popinet. Numerical models for surface tension.
Annu. Rev. Fluid Mech. , 50:49–75, 2018.[67] A. Prosperetti and G. Tryggvason.
Computational Methods for Multiphase Flow . Cambridge UniversityPress, 2007.[68] L. Qian, Y. Wei, and F. Xiao. Coupled thinc and level set method: A conservative interface capturingscheme with high-order surface representations.
J. Comput. Phys. , 373:284–303, 2018.[69] M. Raessi and H. Pitsch. Consistent mass and momentum transport for simulating incompressibleinterfacial flows with large density ratios using the level set method.
Comput. Fluids , 63:70–81, 2012.[70] J.C. Ramirez, C. Beckermann, A.s. Karma, and H.-J. Diepers. Phase-field modeling of binary alloysolidification with coupled heat and solute diffusion.
Physical Review E , 69(5):051607, 2004.[71] F. R¨osler and D. Br¨uggemann. Shell-and-tube type latent heat thermal energy storage: numericalanalysis and comparison with experiments.
Heat and mass transfer , 47(8):1027, 2011.[72] M. Rudman. A volume-tracking method for incompressible multifluid flows with large density varia-tions.
Int. J. Numer. Methods. Fluids , 28:357–378, 1998.[73] M. Salcudean and Z. Abdullah. On the numerical modelling of heat transfer during solidificationprocesses.
International journal for numerical methods in engineering , 25(2):445–473, 1988.[74] Zaki Saptari Saldi.
Marangoni driven free surface flows in liquid weld pools . PhD thesis, Delft Universityof Technology, 2012.[75] A.A. Samarskii, P.N. Vabishchevich, O.P. Iliev, and A.G. Churbanov. Numerical simulation of con-vection/diffusion phase change problems—a review.
International journal of heat and mass transfer ,36(17):4095–4106, 1993. 2876] R. Scardovelli and S. Zaleski. Direct numerical simulation of free-surface and interfacial flow.
Annu.Rev. Fluid Mech. , 31:567–603, 1999.[77] J.A. Sethian and P. Smereka. Level set method for fluid interfaces.
Annu. Rev. Fluid Mech. , 35:341–372,2003.[78] V. Shatikian, G. Ziskind, and R. Letan. Numerical investigation of a pcm-based heat sink with internalfins.
International journal of heat and mass transfer , 48(17):3689–3706, 2005.[79] J Shen. Modeling and numerical approximation of two-phase incompressible flows by a phase-fieldapproach.
Multiscale Modeling and Analysis for Materials Simulation , 22:147–195, 2011.[80] H. Shmueli, G. Ziskind, and R. Letan. Melting in a vertical cylindrical tube: Numerical investigationand comparison with experiments.
International Journal of Heat and Mass Transfer , 53(19-20):4082–4091, 2010.[81] K.R. Sultana, S.R. Dehghani, K. Pope, and Y.S. Muzychka. Numerical techniques for solving so-lidification and melting phase change problems.
Numerical Heat Transfer, Part B: Fundamentals ,73(3):129–145, 2018.[82] M. Sussman, P. Smereka, and S. Osher. A level set approach for computing solutions to incompressibletwo-phase flow.
J. Comput. Phys. , 114:146–159, 1994.[83] W. Tan, N.S. Bailey, and Y.C. Shin. A novel integrated model combining cellular automata and phasefield methods for microstructure evolution during solidification of multi-component and multi-phasealloys.
Computational Materials Science , 50(9):2573–2585, 2011.[84] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop,D. Lifka, G.D. Peterson, R. Roskies, J.R. Scott, and N. Wilkins-Diehr. Xsede: accelerating scientificdiscovery.
Comput. Sci. Eng. , 16:62–74, 2014.[85] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, and Y.J.Jan. A front-tracking method for the computations of multiphase flow.
J. Comput. Phys. , 169:708–759,2001.[86] G. Tryggvason, R. Scardovelli, and S. Zaleski.
Direct Numerical Simulations of Gas-Liquid MultiphaseFlows . Cambridge University Press, 2011.[87] S.O. Unverdi and G. Tryggvason. A front-tracking method for viscous, incompressible, multi-fluidflows.
J. Comput. Phys. , 100:25–37, 1992.[88] J. Vogel and A. Thess. Validation of a numerical model with a benchmark experiment for melt-ing governed by natural convection in latent thermal energy storage.
Applied Thermal Engineering ,148:147–159, 2019.[89] V.R. Voller, M. Cross, and N.C. Markatos. An enthalpy method for convection/diffusion phase change.
International journal for numerical methods in engineering , 24(1):271–284, 1987.[90] V.R. Voller and C. Prakash. A fixed grid numerical modelling methodology for convection-diffusionmushy region phase-change problems.
International Journal of Heat and Mass Transfer , 30(8):1709–1719, 1987.[91] V.R. Voller and C.R. Swaminathan. Eral source-based method for solidification phase change.
Nu-merical Heat Transfer, Part B Fundamentals , 19(2):175–189, 1991.[92] W. Voller. An overview of numerical methods for solving phase change problems.
Advances in numericalheat transfer , 1:341, 1996.[93] F Xiao, Y Honma, and T Kono. A simple algebraic interface capturing scheme using hyperbolic tangentfunction.
Int. J. Numer. Meth. Fluids , 48(9):1023–1040, 2005.2994] B. Xie, Jin P., Du. Y., and S. Liao. A consistent and balanced-force model for incompressible multiphaseflows on polyhedral unstructured grids.
International Journal of Multiphase Flow , 122:103125, 2020.[95] B. Xie and F. Xiao. Toward efficient and accurate interface capturing on arbitrary hybrid unstructuredgrids: The thinc method with quadratic surface representation and gaussian quadrature.
J. Comput.Phys. , 349:415–440, 2017.[96] J. Yan, W. Yan, S. Lin, and G.J. Wagner. A fully coupled finite element formulation for liquid–solid–gas thermo-fluid flow with melting and solidification.
Computer Methods in Applied Mechanics andEngineering , 336:444–470, 2018.[97] W. Yan, W. Ge, Y. Qian, S. Lin, B. Zhou, W.K. Liu, F. Lin, and G.J. Wagner. Multi-physics mod-eling of single/multiple-track defect mechanisms in electron beam selective melting.
Acta Materialia ,134:324–333, 2017.[98] P. Yue, J.J. Feng, C. Liu, and J. Shen. A diffuse-interface method for simulating two-phase flows ofcomplex fluids.
J. Fluid Mech. , 515:293–317, 2004.[99] A. Zhang, J. Du, X. Zhang, Z. Guo, Q. Wang, and S. Xiong. Phase-field modeling of microstructureevolution in the presence of bubble during solidification.
Metallurgical and Materials Transactions A ,51(3):1023–1037, 2020.[100] J. Zhang and X. Yang. A fully decoupled, linear and unconditionally energy stable numerical scheme fora melt-convective phase-field dendritic solidification model.
Computer Methods in Applied Mechanicsand Engineering , 363:112779, 2020.[101] H. Zhao, W. Niu, B. Zhang, Y. Lei, M. Kodama, and T. Ishide. Modelling of keyhole dynamicsand porosity formation considering the adaptive keyhole shape and three-phase coupling during deep-penetration laser welding.