A consistent and conservative model and its scheme for N-phase-M-component incompressible flows
AA consistent and conservative model and its scheme for N -phase- M -component incompressible flows 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
Abstract
In the present work, we propose a consistent and conservative model for multiphase and multicom-ponent incompressible flows, where there can be arbitrary numbers of phases and components. Eachphase has a background fluid called the pure phase, each pair of phases is immiscible, and componentsare dissolvable in some specific phases. The model is developed based on the multiphase Phase-Fieldmodel including the contact angle boundary condition, the diffuse domain approach, and the analyseson the proposed consistency conditions for multiphase and multicomponent flows. The model conservesthe mass of individual pure phases, the amount of each component in its dissolvable region, and thus themass of the fluid mixture, and the momentum of the flow. It ensures that no fictitious phases or com-ponents can be generated and that the summation of the volume fractions from the Phase-Field modelis unity everywhere so that there is no local void or overfilling. It satisfies a physical energy law andit is Galilean invariant. A corresponding numerical scheme is developed for the proposed model, whoseformal accuracy is 2nd-order in both time and space. It is shown to be consistent and conservative andits solution is demonstrated to preserve the Galilean invariance and energy law. Numerical tests indicatethat the proposed model and scheme are effective and robust to study various challenging multiphaseand multicomponent flows.
Keywords:
Multi-phase; Multi-component; Phase-Field model; Diffuse domain approach; Consistentscheme; Conservative scheme
Multiphase and multicomponent flow problems are ubiquitous and have various applications. The problemsinclude strong interactions among phases and components, leading to challenges in developing analytical ornumerical solutions. Before we summarize the related studies, we first clarify the terminology in the presentwork, since we notice that “multiphase” and “multicomponent” are commutable in some literature. In thepresent work, we consider that phases are immiscible with each other and the domain is occupied by at leastone of the phases, and that components are materials dissolvable in some specific phases and their presencecan increase the local density and viscosity based on their amounts (or concentrations).The multiphase problems are more difficult to be modeled numerically since a successful numerical modelshould at least be able to overcome the numerical diffusion and dispersion due to the discontinuity at themoving and deforming phase interfaces. Many efforts have been focused on the two-phase problems inrecent decades and the one-fluid formulation [64, 45] is the most popular framework where the motionsof different phases are modeled by a single momentum equation. Under this framework, the locations ofthe two phases can be specified by various methods, e.g., the Front-Tracking method [65, 63], the Level-Set method [40, 58, 52, 17, 18], the conservative Level-Set method [38, 39, 7], the Volume-of-Fluid (VOF) ∗ 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 ] M a r ethod [20, 49, 43, 41], the ”THINC” method [67, 26, 68, 46], and the Phase-Field (or Diffuse-Interface)method [1, 27, 53, 24]. The effect of the surface tension can be modeled by, e.g., the continuous surfaceforce model (CSF) [6] and the ghost fluid method (GFM) [13, 34], and they are integrated into the flowdynamics by the balanced-force algorithm [16, 44]. Some recent studies start considering the problemsincluding more than two phases and extending the two-phase models to the three- or multi-phase cases with,e.g., the Volume-of-Fluid (VOF) and Moment-of-Fluid (MOF) methods [50, 51, 15, 42], and the Level-Setmethod [55, 37, 57]. Nevertheless, a growing number of studies use a Phase-Field model due to its simplicityand effectiveness for three-phase problems, e.g., [3, 4, 33, 29, 71, 70], and for N -phase problems, e.g.,[5, 30, 35, 31, 66, 32]. Some encouraging progresses on Phase-Field modeling for N -phase flows have beendone by Dong [9, 10, 11, 12], and both the Phase-Field equation and its boundary condition are studied. Thelatest model in [12] is the first Phase-Field model that is fully reduction consistent. Based on the Phase-Fieldmodel in [12] and the consistency analysis [24], Huang et al. developed a consistent and conservative schemefor multiphase incompressible flows [21]. More recently, Huang et al. develop a consistent and conservativevolume distribution algorithm for multiphase flows and apply it to both the Cahn-Hilliard and conservativeAllen-Cahn models [23], so that the models and their numerical solutions are consistent, conservative, andbounded.Despite the active studies on the modeling and simulation of two-/multi-phase flows, the studies includingmultiple components in a multiphase system are relatively rare, probably because the multiphase problemitself is adequately complicated, and a general and physically plausible multiphase model, e.g., the one in[12, 21], is developed very recently. When there are multiple phases and components in a system, a componentcan dissolve in several different phases, and, at the same time, there can be different components present ina single phase. It is challenging to incorporate all these relationships between phases and components in ageneral model. In addition, the phases are moving, deforming, and even experiencing topological changes,and the components inside the phases need to respond to these motions appropriately, which casts anotherchallenge in modeling. The appearance of the components may change the local density, which influences themass transport and the consistency of the model. We notice that some existing works can be categorized intothe two-phase-one-component problems, e.g., the problems including a surfactant [61, 54, 56, 72], where theconcentration of the surfactant changes the surface tension and, as a result, introduces the Marangoni effect,and the problems including phase change [19, 48], where the vapor generated from the liquid is dissolvable inthe condensible gas. However, including these additional physics in the problems having multiple phases andcomponents is out of the scope of the present study. We confine our work to incompressible flows withoutphase change and assume that the Marangoni effect due to the presence of the components is negligible.In the present work, we develop a consistent and conservative model for multiphase and multicomponentincompressible flows. The model allows arbitrary numbers of phases and components appearing simulta-neously. Each component can exist in different phases. In each phase, there is a background fluid calledthe pure phase, and multiple components can be dissolved in this pure phase. Therefore, each phase is a“solution” of its pure phase (or background fluid) as the “solvent” and the components dissolved in thisphase as the “solutes”. Individual pure phases and components have their own densities and viscosities.Each pair of phases has a surface tension and each component has a diffusivity in a given phase. At awall boundary, each pair of phases has a contact angle. The model is based on the multiphase Phase-Fieldmodel including the contact angle boundary condition [11, 12, 21], the diffuse domain approach [36], andthe consistency analysis [24]. The Phase-Field model introduces a set of order-parameters to indicate thelocation of each phase. The interfaces between phases are replaced by interfacial regions whose thickness ispreserved by the balance of the thermodynamical compression and diffusion. The diffuse domain approachis applied to replace the convection-diffusion equation of a component defined in a specific phase with itsmathematical equivalent equation defined in the whole domain. The consistency analysis is performed basedon the consistency conditions for multiphase and multicomponent flows, proposed in the present work. It en-sures that the flow dynamics is correctly modeled for different phases and components, and that no fictitiousphases or components are allowed to be generated. The resulting model conserves the mass of individualpure phases, the amount of each component in its dissolvable region, and thus the mass of the fluid mixture,and the momentum of the multiphase and multicomponent flows. It ensures that the summation of thevolume fractions from the Phase-Field model is unity everywhere so that there is no local void or overfilling.It satisfies a physical energy law and it is Galilean invariant. It satisfies all the consistency conditions, whichare the consistency of reduction, the consistency of volume fraction conservation, the consistency of mass2onservation, and the consistency of mass and momentum transport. It should be noted that the consistencyconditions play a critical role in avoiding fictitious phases and components, in deriving the energy law, andin proving the Galilean invariance of the model. The model is flexible to allow a component to cross a phaseinterface with either zero flux or continuous flux based on the property of the phase interface if the compo-nent is dissolvable in both sides of the interface. In addition to study the dynamics of the multiphase andmulticomponent flows, the present model is applicable to study some multiphase flows, where the miscibilityof each two phases can be different.The corresponding numerical scheme is developed for the proposed model, which is an extension of ourprevious schemes for multiphase flows [21, 23]. The scheme is formally 2nd-order accurate in both timeand space and it preserves the properties of the model. The scheme conserves the mass of individual purephases, the amount of each component in its dissolvable region, and therefore the mass of the fluid mixture.The momentum is exactly conserved by the scheme with the conservative method for the interfacial force,while it is essentially conserved with the balanced-force method. All the consistency conditions are shownto be satisfied in the discrete level by the scheme. This is very important for the problems including largedensity ratios [24] and to ensure that no fictitious phases or components can be generated by the scheme.The numeral solution appears to preserve the Galilean invariance and satisfy the energy law. The propertiesand capabilities of the proposed model and scheme are demonstrated numerically.The rest of the paper is presented as follows. In Section 2, the problem of interest is defined in detail,and the consistent and conservative N -phase- M -component model is introduced, followed by the consistencyanalysis, derivation of the energy law, and the proof of Galilean invariance. In section 3, the numericalscheme for the proposed model is introduced, followed by the analysis and discussion about its accuracy,consistency, and conservation properties. In Section 4, multiple numerical cases are performed to demonstratethe properties of the proposed model and scheme, and to show their capability for handling challengingproblems. Finally, we conclude the present study in Section 5. N -phase- M -component model In this section, we first define the problem of interest. Then, the multiphase Phase-Field model, the dif-fuse domain approach for individual components, the material properties, and the momentum equation areintroduced. The consistency conditions for the N -phase- M -component system are proposed and the consis-tency analysis is performed. Thanks to satisfying the consistency conditions, the proposed model satisfiesa physical energy law and is Galilean invariant. At the end of the section, we summarize the propertiesof the N -phase- M -component model, and provide an example of a multiphase problem where each pair ofphases can be either miscible or immiscible. We then illustrate the applicability of the model in differentcircumstances of cross-interface transports of a component that is dissolvable in both sides of the interface. The problem considered in the present work includes multiple phases and components, and we use N ( N (cid:62) M ( M (cid:62)
0) for the number of components.The N phases are specified by defining a set of order parameters { φ p } Np =1 ∈ [ − , φ p = 1, there is only Phase p , while φ p = − p . The volume fractions of individual phases are related to { φ p } Np =1 by χ p = 1 + φ p , (cid:54) p (cid:54) N, (1)so that { χ p } Np =1 ∈ [0 , N (cid:88) p =1 χ p = 1 , (2)3r equivalently N (cid:88) p =1 φ p = 2 − N, (3)everywhere. Each phase has a background fluid called the pure phase, and ρ φp and µ φp denote the density andviscosity, respectively, of pure Phase p . Every pair of phases is immiscible, and therefore there is a pairwisesurface tension between them, denoted by σ p,q for Phases p and q . The pairwise surface tensions build a N × N matrix which is symmetric and has zero diagonal. At a wall boundary, there is a pairwise contactangle θ Wp,q of Phases p and q , measured on the side of Phase p . Again, the pairwise contact angles builda N × N matrix whose entries satisfy θ Wp,q + θ Wq,p = π . As already mentioned, the Marangoni effect of thecomponents is assumed to be negligible, and consequently, the pairwise surface tensions and contact anglesdo not change due to the appearance of the components.The M components are represented by a set of concentrations { C p } Mp =1 . Individual components havetheir densities { ρ Cp } Mp =1 and viscosities { µ Cp } Mp =1 . Each component can be dissolved inside different phases.On the other hand, there can be multiple components inside a specific phase. To indicate the dissolvabilitiesamong the components and the phases, the dissolvability matrix with dimension M × N is defined as I Mp,q = (cid:26) , if Component p is dissolvable in Phase q, , else . (4)The diffusion coefficient of Component p in Phase q is denoted by D p,q . Noted that values of { C p } Mp =1 aremeaningful only inside the corresponding phases they dissolve in. To obtain a meaningful value in the entiredomain of interest, one needs to use, e.g., I Mp,q χ q C p , to represent the amount of Component p in Phase q .As pointed out in [60], there are two commonly used ways to interpret concentration. The first case is toquantify the amount of a dissolved solute of a solution such as the “molar concentration” of salt in a salt watersolution. The second case is to quantify the composition of miscible fluids such as the “volume fraction” ofethanol in a water-ethanol system. Here, the components act like “solutes” while the pure phases behave as“solvents”. Therefore, the concentrations of the components are interpreted as the first case, e.g., the “molarconcentration”. We further assume that inside each phase the “solution” compound of the pure phase (asthe “solvent”) and the components (as the “solutes”) dissolved in that phase is dilute. In other words, ina specific phase, the volume fractions of the components dissolved in it are negligibly small compared tothe volume fraction of the pure phase or the background fluid of that phase. Therefore, the volume of eachphase will not be changed by either the appearance or cross-interface transport of the components. Thisassumption of diluteness is also consistent with neglecting the Marangoni effect of the components. It shouldbe noted that the concentrations in the proposed model can also be interpreted as the second case, e.g., the“volume fraction”, in a subset of problems where each phase has a single diffusion coefficient shared by allthe components dissolved in it and the cross-interface transport of those components is prohibited. Althoughthe assumption of diluteness can be removed in these problems, one has to be cautious about the assumptionof neglecting the Marangoni effect of the components when applying the proposed model to these problems.In the present work, we loosely call { ρ Cp } Mp =1 and { µ Cp } Mp =1 densities and viscosities, respectively, regardlessof which interpretation of the concentrations is employed. The actual meaning of { ρ Cp } Mp =1 depends on theconcrete definition of the concentrations. For example, ρ Cp is the molar mass of Component p when theconcentrations are “molar concentration”, while it is the density of pure Component p minus the densityof the pure phase the component dissolves in when the concentrations mean “volume fraction”. The sameworks for { µ Cp } Mp =1 .The densities of the pure phases and components are all constant, and phase changes are not consideredin the present work, in addition to the aforementioned assumptions. As a result, the velocity of the flow u is divergence-free, i.e., ∇ · u = 0 . (5) The governing equations of the N -phase- M -component model are described in detail, which include themultiphase Phase-Field model, the diffuse domain approach for the components, the material properties,4nd the momentum equation. The dynamics of the volume fraction contrasts { φ p } Np =1 is governed by the following multiphase Phase-Fieldmodel [12, 21] defined in the whole domain of interest Ω, ∂φ p ∂t + ∇ · ( u φ p ) = ∇ · (cid:32) N (cid:88) q =1 M φp,q ∇ ξ q (cid:33) , (cid:54) p (cid:54) N, (6)where M φp,q is the mobility, and ξ p is the chemical potential of Phase p . The homogeneous Neumannboundary condition is employed to the chemical potentials at the boundary of Ω unless otherwise specified.The definitions of M φp,q and ξ p are M φp,q = M φ ( φ p , φ q ) = (cid:26) − M φ (1 + φ p )(1 + φ q ) , p (cid:54) = q,M φ (1 + φ p )(1 − φ q ) , p = q, (cid:54) p, q (cid:54) N, (7)and ξ p = N (cid:88) q =1 λ p,q (cid:18) η ( g (cid:48) ( φ p ) − g (cid:48) ( φ p + φ q )) + ∇ φ q (cid:19) , (cid:54) p (cid:54) N, (8)along with the mixing energy density λ p,q = 32 √ ησ p,q , (cid:54) p, q (cid:54) N, (9)and with the two potential functions g ( φ ) = 14 (1 − φ ) , g ( φ ) = 14 φ ( φ + 2) , (10)where M φ is a positive constant, η represents the thickness of the interfaces, and g (cid:48) ( φ ) and g (cid:48) ( φ ) are thederivatives of g ( φ ) and g ( φ ) with respect to φ , respectively.Equivalently, the Phase-Field equation, Eq.(6), can be written as ∂φ p ∂t + ∇ · m φ p = 0 , (cid:54) p (cid:54) N, (11)where m φ p = u φ p − N (cid:88) q =1 M φp,q ∇ ξ q , (cid:54) p (cid:54) N, (12)is the Phase-Field flux.At a wall boundary, the contact angle boundary condition in [11, 21], i.e., n · ∇ φ p = N (cid:88) q =1 ζ p,q φ p φ q , (cid:54) p (cid:54) N, (13)where ζ p,q = 2 √ η cos( θ Wp,q ) , (cid:54) p, q (cid:54) N, (14)is implemented, and the contact angles { θ Wp,q } Np,q =1 are set to be π unless otherwise specified in the presentstudy. 5he Phase-Field model Eq.(6) is first developed by Dong [12] in terms of volume fractions and is laterreformulated in its conservative form in terms of volume fraction contrasts by Huang et al. in [21]. Oneimportant property of the Phase-Field model is that it satisfies the consistency of reduction, which is definedin Section 2.3. This property grantees that no fictitious phase is allowed by the model to appear. Each φ p isgoverned by a conservative law, and thus, with a proper boundary condition, e.g., that the normal velocityvanishes at the domain boundary, ddt (cid:82) Ω φ p d Ω = 0 is true for all p . This is equivalent to ddt (cid:82) Ω χ p d Ω = 0 forall p from Eq.(1), which represents the conservation of the phase volumes, and further to ddt (cid:82) Ω ρ φp χ p d Ω = 0for all p , which is the mass conservation of the pure phases. Moreover, it can be shown that the summationconstraint, i.e., Eq.(3), is implied by Eq.(6). Summing Eq.(6) over p , both sides of the summed equationbecome zero, thanks to the definition of M φp,q in Eq.(7), if Eq.(3) is true. In other words, Eq.(3) is the solutionto the summed equation. More details about the Phase-Field model are referred to [12]. Both studies in[12, 21] indicate that the Phase-Field model is effective to model multiphase flows when it is appropriatelycombined with the momentum equation. An alternative Phase-Field model is the conservative Allen-Cahnmodel for multiphase flows developed in [23], which shares the same properties as the model introduced inthe present work. However, it doesn’t include the effect of the contact angles. Suppose Component p is dissolvable in Phase q , i.e., I Mp,q = 1, the dynamics of Component p is modeled bythe convection-diffusion equation, i.e., ∂C p ∂t + ∇ · ( u C p ) = ∇ · ( D p,q ∇ C p ) , (15)defined in Ω q , which is the part of the domain occupied by Phase q . The boundary condition of Component p at Γ q,r , i.e., the boundary between Ω q and Ω r , is n q,r · ( D p,q ∇ C p ) = (cid:26) , if I Mp,r = 0 − n r,q · ( D p,r ∇ C p ) , if I Mp,r = 1 , , q (cid:54) = r, (16)where n q,r is the normal vector pointing from Ω q to Ω r , and it is obvious that n q,r = − n r,q . The boundarycondition Eq.(16) tells that if Component p is unable to dissolve in Ω r , i.e., I Mp,r = 0, no diffusive flux isallowed across Γ q,r . On the other hand, if Component p is also dissolvable in Ω r , i.e., I Mp,r = 1, the diffusiveflux is continuous at Γ q,r . The convective flux is zero at Γ q,r by considering that the velocities of the flowand of Γ q,r are the same. As already mentioned in Section 2.1, the effect of cross-interface transports ofcomponents, i.e., Eq.(16), on changing phase volumes is not counted, thanks to the assumption of diluteness.In addition, by an appropriate setup, the proposed model allows to remove the cross-interface transports ofcomponents without modifying Eq.(16). Related details are in Section 2.6.Since Phase q can be moving, deforming, or even experiencing topological change, Ω q is time-dependentand its boundary can be extremely complex, which casts great challenges on solving Eq.(15) defined in itand on assigning Eq.(16) at its boundary. It is desirable to extend the definition of Eq.(15) from in Ω q to inthe whole domain of interest Ω, which is fixed, and to applied the boundary condition Eq.(16) implicitly likemodeling the surface tension with the continuum surface force. The mathematical equivalence of differentdifferential operators defined in a time-dependent and complex domain to those defined in a larger and fixeddomain is proposed by [36], where some additional terms appear to take the boundary condition into account.Consequently, Eq.(15) defined in Ω q , along with Eq.(16) at { Γ q,r } Nr =1 ,r (cid:54) = q , has the following mathematicallyequivalent form, ∂ ( H q C p ) ∂t + ∇ · ( u H q C p ) = ∇ · ( H q D p,q ∇ C p ) + N (cid:88) r =1 ,r (cid:54) = q n q,r · ( D p,q ∇ C p ) δ q,r , (17)defined in the whole domain of interest Ω, where H q and δ q,r are the indicator functions of Ω q and Γ q,r suchthat, for a smooth function f , (cid:90) Ω q f d Ω = (cid:90) Ω H q f d Ω , (cid:90) Γ q,r f ds = (cid:90) Ω δ q,r f d Ω . (18)6e again have used the equality of the velocity of the flow to that of Γ q,r . The last term in Eq.(17)represents the contribution from the boundary condition at { Γ q,r } Nr =1 ,r (cid:54) = q . A complete derivation fromEq.(15) to Eq.(17) is available in [36] and interested readers should refer to that. Next, we sum Eq.(17) overall the q that I Mp,q = 1, and the resultant equation for Component p is ∂ ( H Mp C p ) ∂t + ∇ · ( u H Mp C p ) = ∇ · ( D p ∇ C p ) , (19)where H Mp = N (cid:88) q =1 ,I Mp,q =1 H q = N (cid:88) q =1 I Mp,q H q , (20) D p = N (cid:88) q =1 ,I Mp,q =1 H q D p,q = N (cid:88) q =1 I Mp,q H q D p,q . (21)Due to the boundary condition Eq.(16), the summation of the last term in Eq.(17) over all the componentsis zero, given that δ q,r = δ r,q . H Mp indicates the dissolvable region of Component p . D p can be understoodas an effective diffusion coefficient of Component p in the whole domain Ω, whose value is zero in Ω r | I Mp,r =0 and D p,q in Ω q | I Mp,q =1 .As shown in [36], the exact indicator function H can be approximated by a smoothed indicator function H ε , which has a transition from 0 to 1 with a thickness ε , and the resulting equation using the smoothedindicator function asymptotically converges to the original equation and boundary condition, e.g., Eq.(15)and Eq.(16) in the present case, as ε tends to zero. Such a method is called the diffuse domain approach,which is developed by Li et al. [36] for solving a partial differential equation (PDE) in a complex domain.As a result, the volume fractions { χ p } Np =1 from the Phase-Field model is a direct candidate of { H εp } Np =1 and ε is the same as η in this case.Before replacing { H p } Np =1 by { χ p } Np =1 , we notice that Eq.(19) implies the dynamic equations for { H p } Np =1 .If C p is homogeneous, Eq.(19) along with Eq.(20) becomes ∂H p ∂t + ∇ · ( u H p ) = 0 , (cid:54) p (cid:54) N. (22)If H p is replaced by χ p , Eq.(22) does not hold, since, from Eq.(1), Eq.(11) and Eq.(12), the dynamics of thevolume fractions is governed by ∂χ p ∂t + ∇ · m χ p = 0 , (cid:54) p (cid:54) N, (23)where m χ p = 12 ( u + m φ p ) , (cid:54) p (cid:54) N, (24)is the volumetric fluxes. As a result, directly replacing { H p } Np =1 by { χ p } Np =1 leads to inconsistency of Eq.(22)with Eq.(23). In order to assure the consistency, the dynamic equations for { χ p } Np =1 , derived from thecomponent equations, should be the same as Eq.(23) instead of Eq.(22), when { χ p } Np =1 from the Phase-Field model is used to approximate { H p } Np =1 . We call such a consistency the consistency of volume fractionconservation, which is important for the energy law and the Galilean invariance of the model, as shown inSections 2.4 and 2.5, and has not been noticed and discussed in the original diffuse domain approach [36].Consequently, the consistent component equation is ∂ ( χ Mp C p ) ∂t + ∇ · ( m χ Mp C p ) = ∇ · ( D p ∇ C p ) , (cid:54) p (cid:54) M, (25)7here χ Mp = N (cid:88) q =1 I Mp,q χ q , (cid:54) p (cid:54) M, (26) m χ Mp = N (cid:88) q =1 I Mp,q m χ q , (cid:54) p (cid:54) M, (27) D p = N (cid:88) q =1 I Mp,q χ q D p,q , (cid:54) p (cid:54) M. (28)Unless otherwise specified, the homogeneous Neumann boundary condition is applied to Eq.(25) at theboundary of domain Ω. It is obvious that Eq.(23) is recovered by Eq.(25) when C p is homogeneous. FromEq.(25), the total amount of Component p inside its dissolvable region χ Mp is conserved, i.e., ddt (cid:82) Ω ( χ Mp C p ) d Ω =0. After defining the component flux as m C p = m χ Mp C p − D p ∇ C p , (cid:54) p (cid:54) M, (29)Eq.(25) can be written as ∂ ( χ Mp C p ) ∂t + ∇ · m C p = 0 , (cid:54) p (cid:54) M. (30) Phase p is composed of its pure phase and the components dissolved in it, both of which contribute to themass in Ω p . As a result, the density of Phase p in Ω p , denoted by ρ p , is ρ p = ρ φp + M (cid:88) q =1 ,I Mq,p =1 ρ Cq C q = ρ φp + M (cid:88) q =1 I Mq,p ρ Cq C q . (31)The density of the fluid mixture in Ω, denoted by ρ , is the volume average of { ρ p } Np =1 , i.e., ρ = N (cid:88) p =1 ρ p χ p = N (cid:88) p =1 ρ φp χ p + M (cid:88) p =1 ρ Cp χ Mp C p , (32)where the first term in the rightmost is contributed from the pure phases (or the background fluids) and thesecond term appears due to the components. Similarly, the viscosity of the fluid mixture in Ω is µ = N (cid:88) p =1 µ φp χ p + M (cid:88) p =1 µ Cp χ Mp C p . (33) The velocity of the flow is governed by the momentum equation ∂ ( ρ u ) ∂t + ∇ · ( m ⊗ u ) = −∇ P + ∇ · [ µ ( ∇ u + ∇ u T )] + ρ g + f s , (34)where m is the consistent mass flux, P is the pressure that enforces the divergence-free constraint Eq.(5), g is the gravity, and f s is the surface force that models the multiphase interfacial tensions.8he consistent mass flux reads m = N (cid:88) p =1 ρ φp m χ p + M (cid:88) p =1 ρ Cp m C p . (35)Applying the consistent mass flux in the inertial term of the momentum equation is critical to obtain theenergy law and the Galilean invariance.From [21], the surface force reads f s = 12 N (cid:88) p =1 ξ p ∇ φ p . (36)It is shown in [21] that Eq.(36) is equivalent to (cid:80) Np,q =1 ∇ · ( λ p,q ∇ φ p ⊗ ∇ φ q ). Thus, the momentum of thefluid mixture is conserved by the momentum equation Eq.(34) if the gravity is neglected. The consistency analysis is of great importance for building up a physical multiphase model, as shown inour previous studies of Phase-Field modeling for two- and multi-phase flows [24, 21, 22]. Specifically, theconsistency analysis in the present study is based on the following consistency conditions for N -Phase- M -Component flows. The consistency conditions are • Consistency of reduction: A N -phase- M -component system should be able to recover the corresponding N (cid:48) -phase- M (cid:48) -component system (1 (cid:54) N (cid:48) (cid:54) N −
1, 0 (cid:54) M (cid:48) (cid:54) M − N − N (cid:48) ) phases and( M − M (cid:48) ) components are absent. • Consistency of volume fraction conservation:
The conservation equation of the volume fractions derivedfrom the component equation should be consistent with the one derived from the Phase-Field equationwhen the component is homogeneous. • Consistency of mass conservation:
The mass conservation equation should be consistent with the onederived from the Phase-Field equation, the component equation, and the density of the fluid mixture.The mass flux m in the mass conservation equation should lead to a zero mass source. • Consistency of mass and momentum transport:
The momentum flux in the momentum equation 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.The consistency of reduction is enriched from the N -phase one by adding the component part. It isimportant because it assures that no fictitious phase or component is generated by the model. Any phases orcomponents that are absent at t = 0 will not appear at ∀ t >
0, without considering any sources or injection,if the model satisfies the consistency of reduction. The consistency of reduction of the Phase-Field model,i.e., Eq.(6), and the terms in the density, i.e., Eq.(32), the viscosity, i.e., Eq.(33), and the mass flux, i.e.,Eq.(35), related to { φ p } Np =1 , has been shown in [12, 21]. The momentum equation is reduction consistent aslong as the density ρ , the viscosity µ , and the mass flux m are consistent as well. Thus, our attention is paidon the part related to the components. For a specific p and from Eq.(25), we can obtain C p ≡ , ∀ t >
0, if C p = 0 initially. As a result, the contribution of C p to the density Eq.(32), the viscosity Eq.(33), and themass flux Eq.(35) disappears. Consequently, the M -component system with the absence of Component p reduces to the corresponding ( M −
1) system. We can repeat the procedure ( M − M (cid:48) ) times so that the M -component system reduces to the corresponding M (cid:48) -component system (0 (cid:54) M (cid:48) (cid:54) M − N -phase- M -component model satisfies the consistency of reduction.The consistency of volume fraction conservation is newly proposed in the present work for the diffusedomain approach, and has been considered in Section 2.2.2 when building up the component equationEq.(25). In the original diffuse domain approach, it presumes that the smoothed indicator function has thesame sharp-interface dynamics as the exact one Eq.(22), i.e., the phase interface is advected by the flow9elocity. However, this is not the case when the smoothed indicator function is chosen from the Phase-Field model, where the sharp interface is replaced by the interfacial region whose thickness is small butfinite. There are allowable thermodynamical compression and diffusion, in addition to the flow advection,inside the interfacial region, and this casts a discrepancy between the Phase-Field model and the originaldiffuse domain approach. The consistency of volume fraction conservation aims to remove the discrepancyby incorporating the thermodynamical effects in the Phase-Field model to the component equation. Thisconsistency condition is essential for the model satisfying the energy law and the Galilean invariance, seeSections 2.4 and 2.5.The consistency of mass conservation and the consistency of mass and momentum transport are identicalto those for two- and multi-phase flows [24, 21, 22]. Based on the definitions of the density of the fluidmixture Eq.(32) and the consistent mass flux Eq.(35), it can be shown that the density is governed by ∂ρ∂t + ∇ · m = N (cid:88) p =1 ρ φp (cid:18) ∂χ p ∂t + ∇ · m χ p (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) + M (cid:88) p =1 ρ Cp (cid:32) ∂ ( χ Mp C p ) ∂t + ∇ · m C p (cid:33)(cid:124) (cid:123)(cid:122) (cid:125) = 0 , (37)thanks to Eq.(23) from the Phase-Field equation Eq.(6) and the component equation Eq.(30). Consequently,the consistency of mass conservation is true with the consistent mass flux defined in Eq.(35). Since theconsistent mass flux Eq.(35) has been applied to the inertial term of the momentum equation Eq.(34), theconsistency of mass and momentum transport is satisfied as well. These two consistency conditions arecritical for problems with large density ratios, as indicated in our previous studies [24, 21, 22], since theyhonor the physical coupling between the mass and momentum. Moreover, they are essential for the modelsatisfying the energy law and the Galilean invariance, see Sections 2.4 and 2.5. The proposed N -phase- M -component model satisfies the following physical energy law. We define the freeenergy density, component energy density, and kinetic energy density to be e F = N (cid:88) p,q =1 λ p,q (cid:18) η ( g ( φ p ) + g ( φ q ) − g ( φ p + φ q )) − ∇ φ p · ∇ φ q (cid:19) , (38) e C = M (cid:88) p =1 γ C p χ Mp C p , (39) e K = 12 ρ u · u , (40)respectively, where { γ C p } Mp =1 are positive dimensional constants so that e C has a unit [J / m ] and we choosethem to be unity without loss of generality. The corresponding energies are the integrals of the energydensities over the domain. It should be noted that the chemical potential of Phase p , i.e., Eq.(8), is thefunctional derivative of the free energy with respect to φ p , i.e., ξ p = δE F δφ p = δδφ p (cid:82) Ω e F d Ω.After multiplying the Phase-Field equation Eq.(6) with the chemical potential Eq.(8) and then summingover all the phases p , we obtain the governing equation for the free energy density, i.e.,12 ∂e F ∂t + N (cid:88) p.q =1 λ p,q ∇ · (cid:18) ∂φ p ∂t ∇ φ q + ∂φ q ∂t ∇ φ p (cid:19) + u · N (cid:88) p =1 ξ p ∇ φ p (41)= 12 N (cid:88) p,q =1 ∇ · (cid:0) ξ p M φp,q ∇ ξ q (cid:1) − N (cid:88) p,q =1 M φp,q ∇ ξ p · ∇ ξ q . γ C p C p and then summing over all the components p , we obtain the governing equation for the component energy density, i.e., ∂e C ∂t + M (cid:88) p =1 ∇ · (cid:18) m χ Mp γ C p C p (cid:19) = M (cid:88) p =1 ∇ · ( γ C p D p C p ∇ C p ) − M (cid:88) p =1 γ C p D p ∇ C p · ∇ C p . (42)Eq.(23) has been used in the derivation.After performing the dot product between u and the momentum equation Eq.(34), we obtain the gov-erning equation for the kinetic energy density, i.e., ∂e K ∂t + ∇ · (cid:16) m u · u (cid:17) = −∇ · ( u P ) + ∇ · (cid:0) µ ( ∇ u + ∇ u T ) · u (cid:1) − µ ( ∇ u + ∇ u T ) : ( ∇ u + ∇ u T ) + u · f s . (43)We have neglected the gravity and used Eq.(37) in the derivation.When we combine Eq.(41), Eq.(42) and Eq.(43) together, and assume that all the fluxes vanish at thedomain boundary, we obtain the energy law for the N -phase- M -component system, i.e., dE T dt = ddt (cid:18) E F + E C + E K (cid:19) = ddt (cid:90) Ω (cid:18) e F + e C + e K (cid:19) d Ω (44)= − (cid:90) Ω N (cid:88) p,q =1 M φp,q ∇ ξ p · ∇ ξ q d Ω − (cid:90) Ω M (cid:88) p =1 γ C p D p ∇ C p · ∇ C p d Ω − (cid:90) Ω µ ( ∇ u + ∇ u T ) : ( ∇ u + ∇ u T ) d Ω . Eq.(44) states that the total energy of the N -phase- M -component system, which includes the free energy,the component energy, and the kinetic energy, is not increasing with time without any external input. Theright-hand side (RHS) of Eq.(44) shows three factors to reduce the total energy of the multiphase andmulticomponent system. The first one comes from the thermodynamical non-equilibrium, contributed bythe Phase-Field model. It should be noted that it is effective only in the interfacial region since the mobilityEq.(7) is zero inside the bulk-phase region. The second one results from the inhomogeneity of the componentsin their dissolvable regions, contributed by the component equation. It should be noted that D p is zero where C p is not dissolvable so there is no contribution from those regions to the total energy. The last one is dueto the viscosity of the fluids, contributed from the momentum equation. To include the effect of the contactangles, i.e., the contact angles are other than π , a wall energy functional has to be introduced. However, sofar, it is still an open question whether the contact angle boundary condition Eq.(13) from [11] correspondsto a multiphase wall functional.At the end of this section, we need to emphasize the significance of the consistency conditions in derivingthe energy law. The left-hand side (LHS) of Eq.(42) is derived from C p (cid:32) ∂ ( χ Mp C p ) ∂t + ∇ · ( m χ Mp C p ) (cid:33) = C p (cid:18) χ Mp ∂C p ∂t + m χ Mp · ∇ C p (cid:19) + C p (cid:32) ∂χ Mp ∂t + ∇ · m χ Mp (cid:33)(cid:124) (cid:123)(cid:122) (cid:125) (45)= χ Mp ∂ C p ∂t + m χ Mp · ∇ C p + C p (cid:32) ∂χ Mp ∂t + ∇ · m χ Mp (cid:33)(cid:124) (cid:123)(cid:122) (cid:125) = ∂ ( χ Mp C p ) ∂t + ∇ · (cid:18) m χ Mp C p (cid:19) . Thanks to the consistency of volume fraction conservation, the volume fraction equation, i.e., the terms withunder-brace in Eq.(45), is recovered. Similarly, the derivation of the left-hand side (LHS) of Eq.(43) is u · (cid:18) ∂ ( ρ u ) ∂t + ∇ · ( m ⊗ u ) (cid:19) = ∂e K ∂t + ∇ · (cid:18) m u · u (cid:19) + 12 u · u (cid:18) ∂ρ∂t + ∇ · m (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) = ∂e K ∂t + ∇ · (cid:18) m u · u (cid:19) . (46)The mass conservation equation, i.e., Eq.(37), is recovered in Eq.(46), thanks to the consistency of mass con-servation and the consistency of mass and momentum transport. If the consistency conditions are violated,11he terms with under-brace in Eq.(45) and Eq.(46) are not zero any more. As a result, some additional termswill be added to the total energy of the multiphase and multicomponent system Eq.(44). In other words,even though all the phases are in their thermodynamical equilibrium, all the components are homogeneousin their dissolvable region, and all the fluids are inviscid, the total energy of the multiphase and multicompo-nent system could still change due to the additional terms introduced by the inconsistency. Such a behaviorof the total energy is not physically plausible. Thus, the consistency conditions are playing an essential rolein the physical behavior of the system. The proposed N -phase- M -component model satisfies the Galilean invariance, thanks to the consistencyconditions. In this section, we call the reference frame Frame ( x , t ), and another frame that is moving withrespect to the reference frame with a constant velocity u Frame ( x (cid:48) , t (cid:48) ). All the variables with (cid:48) are measuredin Frame ( x (cid:48) , t (cid:48) ). The Galilean transform is x (cid:48) = x − u t, t (cid:48) = t, u (cid:48) = u − u , f (cid:48) = f, ∇ (cid:48) f (cid:48) = ∇ f, ∇ (cid:48) f (cid:48) = ∇ f, ∂f (cid:48) ∂t (cid:48) = ∂f∂t + u · ∇ f, (47)where f is a scalar function depending on time and space. With the help of the Galilean transform Eq.(47),the governing equations defined in Frame ( x (cid:48) , t (cid:48) ) are ∇ (cid:48) · u (cid:48) = ∇ · ( u − u ) = 0 , (48) ∂φ (cid:48) p ∂t (cid:48) + ∇ (cid:48) · ( u (cid:48) φ (cid:48) p ) − ∇ (cid:48) · (cid:32) N (cid:88) q =1 M (cid:48) φp,q ∇ (cid:48) ξ (cid:48) q (cid:33) = ∂φ p ∂t + u · ∇ φ p + ∇ · (( u − u ) φ p ) − ∇ · (cid:32) N (cid:88) q =1 M φp,q ∇ ξ q (cid:33) = 0 , (49) m φ (cid:48) p = u (cid:48) φ (cid:48) p − N (cid:88) q =1 M (cid:48) p,q ∇ (cid:48) ξ (cid:48) q = ( u − u ) φ p − N (cid:88) q =1 M p,q ∇ ξ q = m φ p − u φ p , (50) m χ (cid:48) p = 12 ( u (cid:48) + m φ (cid:48) p ) = 12 ( u − u + m φ p − u φ p ) = m χ p − u χ p , (51) m χ (cid:48) Mp = N (cid:88) q =1 I Mp,q m χ (cid:48) q = M (cid:88) q =1 I Mp,q m χ q − u M (cid:88) q =1 I Mp,q χ q = m χ Mp − u χ Mp , (52) ∂ ( χ (cid:48) Mp C (cid:48) p ) ∂t (cid:48) + ∇ (cid:48) · ( m χ (cid:48) Mp C (cid:48) p ) − ∇ (cid:48) · ( D (cid:48) p ∇ (cid:48) C (cid:48) p ) = (53) ∂ ( χ Mp C p ) ∂t + u · ∇ ( χ Mp C p ) + ∇ · (cid:16) ( m χ Mp − u χ Mp ) C p (cid:17) − ∇ · ( D p ∇ C p ) = 0 , m C (cid:48) p = m χ (cid:48) Mp C (cid:48) p − D (cid:48) p ∇ (cid:48) C (cid:48) p = ( m χ Mp − u χ Mp ) C p − D p ∇ C p = m C p − u χ Mp C p , (54) m (cid:48) = N (cid:88) p =1 ρ φp m χ (cid:48) p + M (cid:88) p =1 ρ Cp m C (cid:48) p = N (cid:88) p =1 ρ φp ( m χ p − u χ p ) + M (cid:88) p =1 ρ Cp ( m C p − u χ Mp C p ) = m − u ρ, (55) ∂ρ (cid:48) ∂t (cid:48) + ∇ (cid:48) · m (cid:48) = ∂ρ∂t + u · ∇ ρ + ∇ · ( m − u ρ ) = 0 , (56)12 ( ρ (cid:48) u (cid:48) ) ∂t (cid:48) + ∇ (cid:48) · ( m (cid:48) ⊗ u (cid:48) ) + ∇ (cid:48) P (cid:48) − ∇ (cid:48) · [ µ (cid:48) ( ∇ (cid:48) u (cid:48) + ∇ (cid:48) u (cid:48) T )] − ρ (cid:48) g − f (cid:48) s = (57) ∂ ( ρ ( u − u )) ∂t + u · ∇ ( ρ ( u − u )) + ∇ · (( m − ρ u ) ⊗ ( u − u ))+ ∇ P − ∇ · [ µ ( ∇ ( u − u ) + ∇ ( u − u ) T )] − ρ g − f s = ∂ ( ρ u ) ∂t + ∇ · ( m ⊗ u ) − u (cid:18) ∂ρ∂t + ∇ · m (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) + ∇ P − ∇ · [ µ ( ∇ u + ∇ u T )] − ρ g − f s = 0 . The Phase-Field equation Eq.(49), the component equation Eq.(53), the mass conservation equation Eq.(56),the momentum equation Eq.(57), and the divergence-free condition Eq.(48) in Frame ( x (cid:48) , t (cid:48) ) are identical totheir correspondences, i.e., Eq.(6), Eq.(25), Eq.(37), Eq.(34) and Eq.(48), in Frame ( x , t ). Thus the Galileaninvariance is satisfied. It should be noted that the consistency of mass conservation and the consistency ofmass and momentum transport play critical roles in the Galilean invariance of the momentum equation, seethe under-braced term in Eq.(57). If either of the consistency conditions is violated, the under-braced termin Eq.(57) is non-zero, and consequently the Galilean invariance of the momentum equation is failed.The Galilean invariance implies that a homogeneous velocity field is an admissible solution for the momen-tum equation. Consider the case where u (cid:48) = , the divergence-free condition is satisfied and the momentumequation in Frame ( x (cid:48) , t (cid:48) ) requires that −∇ (cid:48) P (cid:48) + ρ (cid:48) g + f (cid:48) s = , (58)which corresponds to the mechanical equilibrium. From the Galilean transform, we have the same mechanicalequilibrium, i.e., −∇ P + ρ g + f s = , (59)and u = u in Frame ( x , t ). So the divergence-free condition is true in Frame ( x , t ). The momentum equationin Frame ( x , t ) becomes ∂ ( ρ u ) ∂t + ∇ · ( m ⊗ u ) = ∇ · [ µ ( ∇ u + ∇ u T )] (60)The right-hand side (RHS) is zero due to ∇ u = . The left-hand side(LHS) is again zero since theconsistency of mass conservation and the consistency of mass and momentum transport recover the massconservation equation Eq.(37). Consequently, the homogeneous velocity u is an admissible solution. A consistent and conservative N -phase- M -component flow model is proposed, which allows different densitiesand viscosities of individual pure phases and components, different surface tensions between every two phases,and different diffusivities between phases and components. The proposed model satisfies all the consistencyconditions, conserves the mass of individual pure phases, conserves the amount of individual componentsinside their dissolvable regions, and conserves the momentum of the multiphase and multicomponent system.It also satisfies the energy law and the Galilean invariance. These properties of the model are not dependenton the material properties of the phases and components considered.The proposed model can be applied to study some multiphase flows where the miscibilities of each pair ofphases are different. For example, to model a three-phase system where Phases 01 and 02 are miscible whilePhases 01 and 03 are immiscible, as well as Phases 02 and 03, we can define a 2-phase-1-component system,where Phase 1 without Component 1 represents Phase 01, Phase 1 with a specific amount, e.g., 1mol / m , ofComponent 1 represents Phase 02, and Phase 2 represents Phase 03. We can set I M , = 1 and I M , = 0. Asa result, Phase 1, with/without Component 1, is immiscible with Phase 2, and the immisciblities betweenPhases 01 and 03 and between Phases 02 and 03 are properly modeled. Component 1 is only dissolvablein Phase 1. When Phase 1 with Component 1 moves to locations where Phase 1 has no Component 1,Component 1 starts to be transported by both convection and diffusion to those locations, which modelsthe miscible behavior of Phases 01 and 02. However, the proposed model is not ready for the problems13here the miscible phases have significantly different surface tensions with respect to a phase that they areimmiscible with or where a phase is miscible with several other phases that are immiscible with each other,since neither the Marangoni effect due to the presence of the components nor the change of phase volumesdue to the transport of components has been considered.The proposed model is flexible to model the cross-interface transport of a component which is dissolvablein both sides of the interface. Although we consider that the diffusive flux is continuous across the phaseinterfaces if the component is dissolvable in both of the phases, the proposed model works in the scenariowhere the component is not allowed to cross the phase interfaces even though it is dissolvable in both sidesof the interface. For example, Component 01 is dissolvable in both Phases p and q but is unable to crossthe phase interfaces between Phases p and q . We can set Component 1 with I M ,p = 1 and I M ,q = 0, andComponent 2 with I M ,p = 0 and I M ,q = 1. Consequently, Component 1 is only allowed to be dissolved inPhase p but not in Phase q , and Component 2 is the opposite. Neither Component 1 nor 2 can cross thephase interface between Phases p and q . Consequently, the concentration of Component 01 in its dissolvableregion χ M C is represented by χ M C + χ M C .In addition to the N -phase- M -component model, more importantly, the present study proposes a generalframework that physically connects the phases, components, and fluid flows, with the help of the proposedconsistency conditions. Thus, it can be implemented to incorporate many other models for locating the phasesor transporting the components. One can freely choose another physical model, instead of Eq.(6), to governthe order parameters. Also, one can replace the convection-diffusion equation Eq.(15) with other transportequation that more accurately produces the dynamics of the components, and obtain the component equationfollowing the diffuse domain approach and the consistency of volume fraction conservation, as described inSection 2.2.2. These result in changing the definitions of the Phase-Field flux and component flux in Eq.(11)and Eq.(30), while the rest of the equations, e.g., the density of the fluid mixture Eq.(32), the consistentmass flux Eq.(35), and the momentum equation Eq.(34) remain unchanged. The Galilean invariance of themomentum equation and the kinetic energy equation Eq.(43) are always valid, regardless of the concretedefinitions of the Phase-Field flux and component flux which depend on the models for the order parametersand components. However, the free energy, the component energy, and the Galilean invariance of the orderparameters and components will be affected. The numerical scheme for the proposed N -phase- M -component model is introduced in this section. Weuse the numerical operators defined in our previous work [24] for spatial discretizations. Specifically, thecollocated grid arrangement is used, where all the dependent variables are stored at the cell centers, and thefluxes and gradients are stored at the cell faces. In addition, there are normal components of the velocitybeing stored at the cell faces, and the cell-face velocity is divergence-free in the discrete sense using the centraldifference at the cell centers. A variable f at ( x i , y j ) is represented by f i,j or [ f ] i,j , and therefore the volumeof cell ( x i , y j ) is denoted by [∆Ω] i,j . We may skip the subscript if the equations work for every cell. Theconvective operators are approximated by the 5th-order WENO scheme [28], and the Laplace operator and thegradient operator are approximated by the 2nd-order central difference. We denote, for example, the discretegradient operator by ˜ ∇ , to distinguish it from its corresponding continuous one. The linear interpolationfrom the nearest nodal values is represented by ¯ f , and any other interpolations or approximations are denotedby ˜ f . The time derivative is approximated by γ t f n +1 − ˆ f ∆ t , where f n +1 denotes the value of f at time level n + 1, ∆ t is the time step size, and γ t and ˆ f are scheme-dependent. Unless otherwise specified, the 2nd-orderbackward difference scheme is used. Therefore, γ t = 1 . f = 2 f n − . f n − . f ∗ ,n +1 is the extrapolationalong the time direction and the 2nd-order one is f ∗ ,n +1 = 2 f n − f n − . Implementations of different kindsof boundary conditions are also described in detail and we refer interested readers to [24]. For a clearpresentation, we consider only two dimensions. The extension to higher dimensions is straightforward, andall the analyses and discussions in this section are valid for higher-dimensional cases.14 .1 Scheme for the Phase-Field equation The scheme for solving the Phase-Field equation Eq.(6) is elaborated in our previous work [21], whichis modified from the one in [12]. The scheme is formally 2nd-order accurate in both time and space, andsatisfies the consistency of reduction in the discrete level, the mass conservation of individual pure phases, i.e., (cid:80) i,j [ φ np ∆Ω] i,j = (cid:80) i,j [ φ p ∆Ω] i,j for all p , and the summation constraint, i.e., (cid:80) Np =1 φ p = 2 − N , at all thecells and time levels. Special cares have been paid to the convection term in Eq.(6) to avoid generatingfictitious phases, local voids or overfilling. After solving the equation, the consistent and conservativeboundedness mapping algorithm [23] is applied to guarantee that { φ p } Np =1 is in [ − ,
1] while preservingall the aforementioned properties of the solution. The resultant fully-discrete Phase-Field equation is γ t φ n +1 p − ˆ φ p ∆ t + ˜ ∇ · ˜ m φ p = 0 , (cid:54) p (cid:54) N, (61)where ˜ m φ p is the discrete Phase-Field flux of Phase p and its definition is available in [23].Once { φ n +1 p } Np =1 and { ˜ m φ p } Np =1 are obtained from Eq.(61), the volume fractions { χ n +1 p } Np =1 and volu-metric fluxes { ˜ m χ p } Np =1 are computed from Eq.(1) and Eq.(24), respectively, and we can proceed to solvethe component equation Eq.(25). The scheme for solving the component equation Eq.(25) is γ t χ M,n +1 p C n +1 p − (cid:92) χ Mp C p ∆ t + ˜ ∇ · ( ˜ m χ Mp ˜ C ∗ ,n +1 p ) = ˜ ∇ · ( D p ˜ ∇ C n +1 p ) , (cid:54) p (cid:54) M, (62)where { χ Mp } Mp =1 , { ˜ m χ Mp } Mp =1 , and { D p } Mp =1 are computed from Eq.(26), Eq.(27), and Eq.(28), respectively,using { χ n +1 p } Np =1 and { ˜ m χ p } Np =1 from Section 3.1.The 2nd-order backward difference scheme is applied in Eq.(62) for the temporal discretization but C n +1 is replaced by C ∗ ,n +1 in the convection term. Note that the difference between C n +1 and C ∗ ,n +1 is O (∆ t ), sothe scheme is still formally 2nd-order accurate in time. The application of the 2nd-order spatial discretizationresults in formally 2nd-order accuracy of the scheme in space. The scheme Eq.(62) conserves the total amountof Component p in its dissolvable region, i.e., (cid:80) i,j [( χ Mp C p ) n +1 ∆Ω] i,j = (cid:80) i,j [( χ Mp C p ) ∆Ω] i,j for all p , if thereis no source of the components at the boundary of the domain. After multiplying ∆Ω to Eq.(62) and thensumming it over all the cells, both the convection and diffusion terms vanish due to the property of thediscrete divergence operator, see the proof in [21, 24]. As a result, (cid:80) i,j [( χ Mp C p )∆Ω] i,j doesn’t change astime advances. The consistency of reduction in the discrete level is also satisfied by the scheme Eq.(62).Suppose C p is zero up to time level n and there is no source of the component at the boundary of the domain,then the convection term and (cid:92) χ Mp C p vanish. The resulting equation is a discretized homogeneous Helmholtzequation. Thus, C p remains to be zero at time level n + 1. In other words, if Component p is absent at t = 0and it doesn’t have any sources at the boundary of the domain, it won’t appear at ∀ t >
0. Therefore theconsistency of reduction is satisfied.The discrete component fluxes are˜ m C p = ˜ m χ Mp ˜ C ∗ ,n +1 p − D p ˜ ∇ C n +1 p , (cid:54) p (cid:54) M, (63)and it also satisfies the consistency of reduction since it becomes zero if C p is zero up to time level n + 1.Once { χ M,n +1 p C n +1 p } Mp =1 and { ˜ m C p } Mp =1 are available, the density, viscosity and mass flux are computedfrom Eq.(32), Eq.(33) and Eq.(35), respectively, and we can proceed to solve the momentum equation. The scheme for the momentum equation with the divergence-free constraint is the same as those in [24] andit has been successfully applied to two- and multi-phase incompressible immiscible flows [24, 21, 22, 23].It is a projection scheme with formally 2nd-order accuracy in both time and space, where the divergence-free constraint is enforced exactly by the cell-face velocity, i.e., ˜ ∇ · u = 0. The momentum is conserved,15.e., (cid:80) i,j [( ρ u ) n ∆Ω] i,j = (cid:80) i,j [( ρ u ) ∆Ω] i,j , if the interfacial tension is not counted, and the momentumconservation does not depend on the explicit form of the mass flux in the inertial term. The interfacial force,i.e., f s in Eq.(36), can be computed by either the balanced-force method, which achieves a better mechanicalequilibrium in the discrete level, or the conservative method, which fully conserves the momentum equation[21]. The proof of the momentum conservation of the scheme is available in [21]. The scheme satisfiesthe consistency of reduction and recovers exactly the fully-discretized single-phase Navier-Stokes equationinside each pure phase if the discrete mass flux ˜ m is reduction consistent. This is true since both { ˜ m φ p } Np =1 and { ˜ m C p } Mp =1 are reduction consistent and the discrete mass flux is computed from them with Eq.(35).The consistency of mass and momentum transport requires that the mass flux in the inertial term of themomentum equation should be the one satisfying the consistency of mass conservation. This statementshould be held at the discrete level. In the following, we show that the discrete mass flux ˜ m from Eq.(35)using { ˜ m φ p } Np =1 and { ˜ m C p } Mp =1 satisfies the consistency of mass conservation.Based on the definitions of volume fractions Eq.(1) and volumetric fluxes Eq.(24), the fully-discretePhase-Field equation Eq.(61), and the divergence-free constraint, we obtain γ t χ n +1 p − ˆ χ p ∆ t + ˜ ∇ · ˜ m χ p = 12 γ t − γ t (cid:122)(cid:125)(cid:124)(cid:123) ˆ1∆ t (cid:124) (cid:123)(cid:122) (cid:125) + ˜ ∇ · u (cid:124) (cid:123)(cid:122) (cid:125) + γ t φ n +1 p − ˆ φ p ∆ t + ˜ ∇ · ˜ m φ p (cid:124) (cid:123)(cid:122) (cid:125) = 0 , (64)which is the discrete counterpart of Eq.(23). As a result, the consistency of volume fraction conservation issatisfied in the discrete level by noticing that Eq.(64) can be derived from Eq.(62) with homogeneous C p .Based on the definition of the discrete component fluxes Eq.(63) and the fully-discrete component equationEq.(62), we obtain χ M,n +1 p C n +1 p − (cid:92) χ Mp C p ∆ t + ˜ ∇ · ˜ m C p = χ M,n +1 p C n +1 p − (cid:92) χ Mp C p ∆ t + ˜ ∇ · ( ˜ m χ Mp ˜ C ∗ ,n +1 p ) − ˜ ∇ · ( D p ˜ ∇ C n +1 p ) = 0 , (65)which is the discrete counterpart of Eq.(30). Thanks to Eq.(64) and Eq.(65), we obtain γ t ρ n +1 − ˆ ρ ∆ t + ˜ ∇· ˜ m = N (cid:88) p =1 ρ φp (cid:32) γ t χ n +1 p − ˆ χ p ∆ t + ˜ ∇ · ˜ m χ p (cid:33)(cid:124) (cid:123)(cid:122) (cid:125) + M (cid:88) p =1 ρ Cp (cid:32) γ t χ M,n +1 p C n +1 p − (cid:92) χ Mp C p ∆ t + ˜ ∇ · ˜ m C p (cid:33)(cid:124) (cid:123)(cid:122) (cid:125) = 0 , (66)which is the discrete counterpart of the mass conservation equation Eq.(37). Therefore, the consistency ofmass conservation is satisfied by the discrete mass flux ˜ m , and thus, the consistency of mass and momentumtransport is satisfied in the momentum equation as well. As discussed in Section 2.5, a constant vector fieldshould be an admissible solution of the momentum equation if the mechanical equilibrium Eq.(59) is satisfied.This is held, independent of the cell size, the time step size, the numbers of phases and components, and thematerial properties of them, by the present scheme for the momentum equation, as long as the consistencyof mass conservation and the consistency of mass and momentum transport are satisfied in the discrete level.The proof is available in [21]. The proposed scheme is decoupled, semi-implicit, and formally 2nd-order accurate in both time and space. Itconserves the mass of individual pure phases, the total amount of the components inside their correspondingdissolvable regions, and, as a result, the total mass of the fluid mixture. The momentum is conserved if theconservative method for the interfacial forces is applied. All the consistency conditions, i.e., the consistencyof reduction, the consistency of volume fraction conservation, the consistency of mass conservation, and theconsistency of mass and momentum transport, are satisfied in the discrete level. The scheme also eliminatesany generation of local void or overfilling. Therefore, the scheme is consistent and conservative. The Galileaninvariance and the energy law of the discrete system are to be explored numerically in the next section.16ne of the advantages of the decoupled semi-implicit scheme is that one solves the complicated systempart by part and only needs to solve linear systems in every part. This avoids solving a huge non-linearcoupled system from a fully implicit scheme and is probably more efficient in a single time step. However,the scheme is only conditionally stable and the time step size ∆ t is restricted. Since all the convection termsare treated explicitly, the CFL condition needs to be considered, i.e., ∆ t (cid:54) ∆ t CFL ∼ hU m where h is thegrid size and U m is the scale of the maximum velocity. Another restriction comes from the surface tension[6, 16], i.e., ∆ t (cid:54) ∆ t σ ∼ (cid:113) h π min( ρ p + ρ q σ p,q ) p (cid:54) = q . The dominant part of the viscous term, i.e., ∇ · ( µ ∇ u ), hasbeen treated implicitly, which greatly removes the dependency of stability on the Reynolds number. Asstated in [14], explicitly treating the remaining part of the viscous term is reasonable because it contributesin only a minor way to the viscous force, and a more careful analysis is available in [2]. The aforementionedcriteria have been commonly used in computational fluid dynamics, and interested readers can refer to[14, 6, 16, 58, 59, 2, 8, 12]. The effect of M φ in Eq.(7) on the stability of the scheme has been discussed in[12], and it is observed that a large M φ could lead to instability. The choice of M φ in the present study isbased on the practices in [12, 21], so that it won’t be the dominant factor of numerical instability but providesreasonably good results. In addition to stability, accuracy is another important aspect when determiningthe time step size, since most interests are in the dynamical behavior of the multiphase and multicomponentproblems. Therefore, a smaller time step size is necessary to produce accurate dynamical data. The efficiencyof the scheme is not a major concern at the moment. Most attention is paid on developing a scheme thathonors the physical properties of the proposed model. In this section, we first investigate the convergence behaviors of the component equation Eq.(25), derivedfrom the diffuse domain approach [36], and its numerical scheme Eq.(62). The other part of the proposedscheme, i.e., the one excluding the component equation, has been analyzed and numerically validated in ourprevious works [24, 22, 21, 23] and we don’t repeat them here. Then, we demonstrate that the model iscapable of accurately capturing different circumstances of cross-interface transports of a component that isdissolvable in both sides of the interface. With the help of the rising bubble problem in [25], the abilityof the proposed model and scheme to produce sharp-interface solution is demonstrated, the effect of theconsistency of volume fraction conservation is illustrated, and the formal order of accuracy of the entirescheme and its convergence to sharp-interface solution is quantified. The Galilean invariance in the discretelevel is validated by a 4-phase-2-component setup including significant differences of material properties.The mass conservation, momentum conservation, and energy law at the discrete level are discussed with thehorizontal shear layer. Finally, two problems are performed to demonstrate the capability of the proposedmodel and scheme in complicated and challenging problems including various critical factors of multiphaseand multicomponent flows.Unless otherwise specified, we chose η = 0 .
01 and M φ = 10 − , use h to denote the cell/grid size, apply thebalanced-force method to discretize the interfacial force in the momentum equation, set the initial velocityto be zero, and call χ Mp C p , which represents the concentration of Component p in its dissolvable region, theconcentration of Component p for convenience. The formal order of accuracy of the scheme for the multiphase incompressible flows, including the Phase-Field equation Eq.(6) and the momentum equation Eq.(34), and the convergence behavior of its numericalsolution to the sharp-interface one have been analyzed and numerically investigated in our previous works[24, 22, 21, 23]. The scheme is formally 2nd-order accurate in both time and space, and the numericalsolution of the model converges to the sharp-interface solution with an accuracy between 1.5th- to 2nd-order. The present case focuses on the newly introduced component equation Eq.(25) derived from thediffuse domain approach [36]. Li et al. [36] used the matched asymptotic analysis based on the interfacialthickness η and showed that the diffuse domain approach recovers, at the leading order, the original equationdefined in a time-dependent domain and the boundary condition at that domain boundary. Their numericalresults showed the convergence of the approach to the solution of the original problem. In this section, we17igure 1: Results of the convergence tests at t = 0 .
05. a) Concentration profiles. b) Errors for different gridsizes.demonstrate the convergence of the diffuse domain approach using a diffusion problem with an analyticalsolution.The domain considered is [1 ×
1] with periodic boundaries along the x axis and with free-slip boundariesalong the y axis. The number of cells along each axis is ranging from 16 to 256. The number of time stepsperformed is 40 for cell size 1 /
16 and is doubled when the cell size is halved. As a result, the time stepsize is proportional to cell size. Two inviscid phases with matched density are considered. Phase 1 is at thebottom below y = 0 . .
1. Its concentration is zero initially, and at the bottom wall, its concentration is constantlyequal to 1.The concentration of Component 1 follows the diffusion equation with a Dirichlet boundary at y = 0and a homogeneous Neumann boundary at y = y . Such a diffusion problem can be solved by separation ofvariables and the exact solution is C S ( y, t ) = C (cid:32) ∞ (cid:88) r =0 − λ r L exp (cid:0) − Dλ r t (cid:1) sin ( λ r ( y − y )) (cid:33) , λ r = 1 + 2 r L π, (67)where y and y are the locations of the Dirichlet and homogeneous Neumann boundaries, L is the lengthof the domain, i.e., L = y − y , and C is the concentration at the Dirichlet boundary. Specifically in thepresent case, y = 0, y = y = 0 . C = 1. It should be noted that Eq.(67) works in [ y , y ], whichis [0 , . ,
1] is H C S , where H isthe exact indicator function of Phase 1 and it is H ( y − y ) with H ( y ) the Heaviside function. However, inthe diffuse domain approach, the exact indicator function H is replaced by the smoothed indicator function χ M , which is (cid:16) (cid:16) y − y √ η (cid:17)(cid:17) in our case. As a result, the best solution to be expected from the diffusedomain approach is χ M C S , and we call it the semi-sharp-interface solution. In the following, we are going tocompare the numerical results to both the sharp-interface solution (S) and the semi-sharp-interface solution(SS) under different circumstances.In the first case, we consider the moment at t = 0 .
05, at which the edge of the boundary layer of theconcentration is far away from the phase interface. Thus, the phase interface should have no effect onthe solution, and the sharp-interface and semi-sharp-interface solutions are identical, i.e., H C S = χ M C S .Under this circumstance, the accuracy of the numerical solution should be the same as its formal order ofaccuracy, which is 2nd-order accurate in both time and space. Our results, shown in Fig.1 match the analysisvery well. Fig.1 a) shows the concentration profiles of the sharp-interface solution H C S , the semi-sharp-interface solution χ M C S and the numerical solution χ M C , and they overlap with each other. The edge18igure 2: Results in the convergence tests at t = 1. a) Profiles of the concentrations, blue solid line: sharp-interface solution, red solid lines: semi-sharp-interface solutions with different cell sizes, yellow dashed lines:numerical solutions with different cell sizes, black solid line: interface. b) Errors for different grid sizes.of the boundary layer of the concentration is at about y = 0 .
3, which is far enough away from the phaseinterface at y = 0 . b) shows the convergence behavior, quantified by the L ∞ and L errors. The point-wise error is defined as the numerical solution minus the sharp-interface or the semi-sharp-interface solution. The L ∞ error is the maximum of the absolute point-wise error and the L error is theroot-mean-square of the point-wise error. The superscripts “S” and “SS” refer to the sharp-interface andsemi-sharp-interface solutions, respectively. We can observe that both L ∞ and L errors are converging tozero with 2nd-order accuracy, as expected. The differences between L S ∞ and L SS ∞ and between L S and L SS are indistinguishable, which implies that the sharp-interface solution is the same as the semi-sharp-interfacesolution. Since the time step size is propositional to the cell size, we can conclude that the scheme is formally2nd-order accurate in both time and space.In the second case, we consider the moment at t = 1, at which the edge of the boundary layer of theconcentration has reached the phase interface. As a result, the sharp-interface solution and the semi-sharp-interface solution are different. We chose η equal to h , so that the numerical solution converges to thesemi-sharp-interface and to the sharp-interface solutions simultaneously as h tends to zero. Fig.2 a) showsthe profiles of the sharp-interface solution, the semi-sharp-interface solution, and the numerical solutionswith different cell sizes. It is clear that the numerical solution, as well as the semi-sharp-interface solution,is approaching the sharp-interface solution, as h and η are getting smaller. Fig.2 b) shows the norms ofthe errors and we observe that L SS ∞ is converging to zero with 2nd-order accuracy and L SS converges evenfaster with 2.5th-order accuracy. However, L S ∞ doesn’t converge to zero while the convergence rate of L S isvery slow at about 0.5th-order accuracy. Since the numerical solution converges to the semi-sharp-interfacesolution, the difference between the semi-sharp-interface solution and the sharp-interface solution dominantlycontributes to the convergence behaviors of L S ∞ and L S . However, those behaviors of L S ∞ and L S are notconsistent with what is observed in Fig.2 a) . To explain this, we use Fig.3, where the sharp-interfacesolution and the semi-sharp-interface solutions with h = η = 0 .
01 and with h = η = 0 .
001 are shown.It is obvious that the difference between the sharp-interface solution and the semi-sharp-interface solutionmajorly appears at the interfacial region, and the difference is much smaller after h and η are 10 timessmaller. However, when we compute the L ∞ errors between the two solutions under the corresponding cellsizes, they are 0 . h = η = 0 .
01 and 0 . h = η = 0 . b) , the maximum difference between the sharp-interfaceand semi-sharp-interface solutions is at the grid point closest to the interface and its location is denoted by y ( h ) s . Although the point-wise error at y (0 . s with h = η = 0 .
001 is smaller than that with h = η = 0 . y (0 . s is closer to the interface than y (0 . s . Consequently, the point-wise error at y (0 . s is larger than thatat y (0 . s with h = η = 0 . y (0 . s with h = η = 0 .
01. As a result, the L ∞ error19igure 3: The sharp-interface and semi-sharp-interface solutions with h = η = 0 .
01 and h = η = 0 .
001 inthe convergence tests at t = 1. a) Concentration profiles. b) Concentration profiles close to the interface.Figure 4: The point-wise errors between the numerical and sharp-interface solutions at y (1 / s in the con-vergence tests at t = 1.between the sharp-interface and semi-sharp-interface solutions doesn’t change much as h and η reduce. Dueto the behavior observed in Fig.3, we consider the convergence of the point-wise error between the numericaland sharp-interface solutions at a specific location. Fig.4 shows the point-wise error at y (1 / s , and by usingthis measure, the error converges to zero with at least 2nd-order accuracy. Two diffusion problems are performed to demonstrate that the proposed model is capable of modelingdifferent circumstances of cross-interface transports of a component that is dissolvable in both sides of theinterface. The domain considered is [1 ×
1] with periodic boundaries along the x axis and free-slip boundariesalong the y axis. The domain is discretized by 128 ×
128 cells and the time step size is ∆ t = 10 − . Phase 1,whose pure phase density is 1000 and viscosity is 10 − , is at the bottom below y = 0 .
3, while the rest of thedomain is filled by Phase 2, whose pure phase density and viscosity are 1 and 2 × − , respectively. Thesurface tension between the phases is 0 . − is dissolvable in both of the phases with diffusivity 0 .
02 inside Phase 1 and 0 .
08 insidePhase 2. Initially, there is no Component 01 inside the domain, while its concentrations are 1 at the bottomwall and 0 . t = 0 . , . , , , , , , ,
10, black solid line: interface. a) results of the first caseof the two diffusion problems. b) results of the second case of the two diffusion problems.In the first case, Component 01 is unable to cross the phase interface between Phases 1 and 2. Asdiscussed in Section 2.6, we need to define Component 1 and Component 2 such that Component 1 isdissolvable only in Phase 1 and Component 2 is dissolvable only in Phase 2. Consequently, the concentrationof Component 01 is represented by χ M C = χ M C + χ M C . The sharp-interface steady solution is a stepfunction jumping from 1 to 0 . y = 0 .
3. In the second case, Component 01 is able to cross the phaseinterface between Phases 1 and 2. Then we only need to define a single component, i.e., Component 1,sharing the same material properties with Component 01 and dissolving in both Phases 1 and 2. As a result,the concentration of Component 01 is χ M C = χ M C . The sharp-interface steady solution is two straightsegments with different slopes connecting at y = 0 .
3, i.e., C S ( y ) = (cid:26) − y + 1 , y (cid:54) . − ( y −
1) + 0 . , y ≥ . χ M C at selected moments in the two cases are shown in Fig.5, along with the corre-sponding sharp-interface steady solutions. As time goes on, the numerical solutions successfully approachthe corresponding sharp-interface steady solutions, even though the background phases have a density ratioof 1000 and a viscosity ratio of 50. Due to the diffuse domain approach, we can observe that the profileschange smoothly across the interface instead of a sharp transition. The results in this section demonstratethat the proposed model is capable of modeling whether a component is allowed to cross a phase interfacewhen it is dissolvable in both sides of the interface. In this section, the rising bubble problem in [25] is considered to demonstrate the capability of the proposedmodel and scheme of producing sharp-interface solutions, to illustrate the effect of the proposed consistencyof volume fraction conservation, and to quantify the convergence of the entire scheme. The domain consideredis [1 × . , .
5) with a radius 0 .
25, while the other part of the domain is filled with a liquid.The bubble has a density 1 and a viscosity 0 .
1, while the density and viscosity of the liquid are 1000 and10, respectively. Their surface tension is 1 .
96, and the gravity is 0 .
98, pointing downward. The motion ofthe bubble is quantified by the trajectory of its center, i.e., y c = (cid:82) Ω yχ b d Ω / (cid:82) Ω χ b d Ω, and the rising speed,i.e., v c = (cid:82) Ω vχ b d Ω / (cid:82) Ω χ b d Ω, where χ b is the volume fraction of the bubble and v is the velocity along the y direction. Here, the mid-point rule is used to compute the integrals. Additional two components are added.Component 1, having a density 0 . .
05, is only dissolvable in Phase 1 with a diffusivity 0 . .
05 and has a density 1600 and viscosity 14.21igure 6: Results of the rising bubble problem under Setups 1 and 2. a) The trajectory of the bubble centerversus time. b) The rising speed of the bubble versus time. The sharp-interface results are from [25].We first compare the results from the proposed model and scheme with the sharp-interface one in [25]using the following two setups. • Setup 1:
The bubble is represented by Phase 1, while the liquid is Phase 2. The components areabsent, by giving their concentrations to be zeros at t = 0. Therefore, the densities, viscosities, andsurface tension of the phases are the same as the given values for the bubble and liquid, and χ b equals χ . The domain is discretized by a grid size of h = . The time step size is ∆ t = 0 . h and η isthe same as h . • Setup 2:
The bubble is represented by Phase 1 with Component 1 having a homogeneous concentration1 dissolved in it. On the other hand, the liquid is composed of Phase 2 and Component 2 with ahomogeneous concentration 0.5. Thus, the concentrations are initially 1 and 0 . . .
05 for theviscosity, while they are 200 for the density and 3 for the viscosity of Phase 2. The rest is the same asSetup 1.As pointed out in Section 2.6, any pure phase can be modeled as another pure phase dissolving componentswith homogeneous concentrations. For example, one can model a salt water as a pure phase or as a solutionof water, which is the pure phase (or the background fluid), and salt, which is the component, with ahomogeneous concentration. Therefore, Setups 1 and 2 should produce the same results. Fig.6 shows theresults from Setups 1 and 2, along with the sharp-interface results from [25]. The results from both setups areindistinguishable from each other and agree well with the sharp-interface results. In addition, under Setup1, the concentrations of the components, which are zero initially, are always zero during the computation.This confirms the consistency of reduction of the model and scheme, as discussed in Sections 2.3 and 3.2,respectively.Next, the effect of the newly proposed consistency condition, i.e., the consistency of volume fractionconservation is considered. The effects of the other consistency conditions have been discussed in, e.g.,[24, 22, 21, 23] for the consistencies of mass conservation and mass and momentum transport, and [5, 35,11, 12, 47, 21, 23] for the consistency of reduction. Setup 2 is repeated but the concentrations of thecomponents are solved from Eq.(19), instead of the proposed Eq.(25), and { H Mp } Mp =1 are directly replacedby { χ Mp } Mp =1 , like the original diffuse domain approach [36], without considering the consistency of volumefraction conservation. Results from solving Eq.(19) are labeled as “IC” which means inconsistent. SinceSetups 1 and 2 are equivalent, the profile of χ Mp C p should be the same as χ Mp ( C p | t =0 ). Fig.7 shows thedifference between χ M C and χ M at x = 0 . χ M C and χ M under Setup 2 at x = 0 .
5. IC: the concentrations are solvedfrom Eq.(19) with { H Mp = χ Mp } Mp =1 . Proposed: the concentrations are solved from the proposed Eq.(25)that satisfies the consistency of volume fraction conservation.growing with time, which probably becomes an origin of numerical instability in highly dynamical problems.On the other hand, the difference is less than 10 − at t = 1 from the proposed consistent model and scheme.The concentration of Component 2 behaviors in a similar manner, and therefore they are not shown here.Finally, the convergence of the entire scheme is considered using Setups 3 and 4, modified from Setup 2,such that the concentrations of the components are not homogeneous initially. • Setup 3 : The concentration of Component 1 is unity inside a circle at (0 . , .
5) with a radius 0 . . y = 0 .
85 but zero elsewhereinitially. The grid size is ranging from to , and η is fixed to be . The solution from the finestgrid is considered as the reference solution. • Setup 4 : The same as Setup 3 but η is changing with the grid size, i.e., η = h .The solutions from different grid sizes in Setup 3 are approximating the same exact solution of the modelwith a fixed η . Therefore, its convergence behavior represents the formal order of accuracy of the proposedscheme. On the other hand, the solutions in Setup 4 converge to a series of solutions that approach the onewith zero η , i.e., the sharp-interface solution. It usually converges slower than the formal order of accuracyof the scheme but is more related to practical applications. Fig.8 shows the results of Setups 3 and 4, andconvergences of the results are evident after successively refining the grid size under both setups. Fig.9 showsthe L norms, i.e., the root mean square, of y c and v c minus their corresponding reference solutions fromthe finest grid. The convergence under Setup 3 is close to but better than 2nd order, which validates theformally 2nd-order accuracy of the proposed scheme. The convergence in Setup 4 is a little slower than thatin Setup 3 but still close to 2nd order. Similar studies have been performed using the multiphase model, i.e.,the one proposed in the present work without the component part, in [21, 23], and the present results areconsistent with those studies, as well as those in Section 4.1 of the present work.The present case demonstrates that the proposed model and scheme are able to produce sharp-interfacesolutions. The proposed consistency of volume fraction conservation is essential to eliminate unphysicalfluctuations of components around phase interfaces. The formal order of accuracy of the entire scheme is2nd-order, and the convergence behavior towards the sharp interface solution is reasonably satisfactory witha rate close to 2nd order. 23igure 8: Results of the rising bubble problem under Setups 3 and 4 using different grid sizes. a) Thetrajectory of the bubble center versus time under Setup 3. b) The rising speed of the bubble versus timeunder Setup 3. c) The trajectory of the bubble center versus time under Setup 4. d) The rising speed of thebubble versus time under Setup 4. In this section, we perform a numerical experiment to validate that the numerical results from the proposedscheme is Galilean invariant. The domain considered is [1 ×
1] with double-periodic boundaries. The domainis discretized by [100 × t = 10 − . We consider four phases and twocomponents. Phase 1, whose pure phase density is 10000 and viscosity is 10 − , is inside a circle of radius0 .
15 at (0 . , . − , is inside a circle of radius0 . . , . − , is inside a horizontalchannel between 0 . (cid:54) y (cid:54) .
4. Phase 4, whose pure phase density is 1 and viscosity is 10 − , fills the rest ofthe domain. The interfacial tensions are σ , = 0 . σ , = 0 . σ , = 0 . σ , = 0 . σ , = 0 . σ , = 0 . − , isdissolvable in Phases 2 and 4, with diffusivity 10 − and 10 − , respectively. It is initially inside a circle ofradius 0 .
05 at (0 . , .
55) with a homogeneous concentration 1. Component 2, whose density is 1 and viscosityis 10 − , is dissolvable in Phases 3 and 4, with diffusivity 10 − and 10 − , respectively. It is initially inside acircle of radius 0 .
05 at (0 . , .
3) with a homogeneous concentration 1. It should be noted that the densityratio considered in the present case is 10 ,
000 and the viscosity ratio is about 1 , − u = −{ , } ,24igure 9: The L norms of y c and v c minus their corresponding reference solutions form the finest grid versusthe grid sizes under Setups 3 and 4.which is equivalent to give u as the initial velocity in this case. The second observer observes both diffusionsof Components 1 and 2, and translations of all the interfaces without deformations. The configurationsobserved by the second observer is labeled as “V”. With given u , the two configurations “S” and “V”should be identical at t = 1 due to the Galilean invariance.Fig.10 shows the configurations of different phases and components at t = 0 and at t = 1 observed bythe two observers. In configuration “S”, the interfaces neither move nor deform, and the components arediffusing only in their corresponding dissolvable regions. The locations of the interfaces at t = 1 is on topof those at t = 0. Component 1, which is initially in Phase 4 only, has diffused into Phase 2 but not toPhases 1 and 3, since Component 1 is not dissolvable in Phases 1 and 3. Component 2, which is dissolvablein Phases 3 and 4, has diffused across the interface between Phases 3 and 4 but not to Phases 1 and 2. Inconfiguration “V”, the interfaces at t = 1 return to their original locations without any deformation, and thebehaviors of the components are the same as those in configuration “S”. The results in configuration “V”has little difference from those in configuration “S”.Fig.11 and Fig.12 further demonstrate the Galilean invariance of the results. Fig.11 shows the profilesof Phase 1 at x = 0 .
3, Phase 2 at x = 0 .
75, and Phase 3 at x = 0 .
9, at t = 0 and t = 1 from bothconfigurations “S” and “V”. The corresponding profiles are overlapping with each other. Fig.12 shows theprofiles of Component 1 at x = 0 . x = 0 . t = 1 from both configurations “S” and“V”. Again, the corresponding profiles are indistinguishable from each other. Consequently, the Galileaninvariance is preserved by the scheme. We discuss the mass conservation, momentum conservation, and energy law using the result from the horizon-tal shear layer. The domain considered is [1 ×
1] with double-periodic boundaries. The domain is discretizedby [128 × t = 0 . h . Phase 1, whose pure phase density is 10 andviscosity is 10 − , is initially in 0 . (cid:54) y (cid:54) .
75. Phase 2, whose pure phase density is 1 and viscosity is 10 − ,fills the rest of the domain. The surface tension between them is 0 . η is √ in this case. Component 1, whose density is 20 and viscosity is 10 − , is dissolvable in Phase 2 only, withdiffusivity 10 − . It initially occupies 0 (cid:54) y (cid:54) .
125 and 0 . (cid:54) y (cid:54)
1, with a homogeneous concentration1. The horizontal velocity is initially 1 inside Phase 1 while it is − t = 0. Middle row:configuration S at t = 1. Bottom row: configuration V at t = 1. Left column: conflagrations of the 4 phases,blue: Phase 1, green: Phase 2, yellow: Phase 3, and White: Phase 4. Middle column: configuration ofComponent 1. Right column: configuration of Component 2. In the middle and right columns, blue solidline: interface of Phase 1 at t = 0, orange dashed line: interface of Phase 1 at t = 1, green solid line: interfaceof Phase 2 at t = 0, red dashed line: interface of Phase 2 at t = 1, yellow solid line: interface of Phase 3 at t = 0, purple dashed line: interface of Phase 3 at t = 1.balanced-force method (B) and the conservative method (C) to discretize the interfacial force.Fig.13 shows the time histories of the total volumes of individual phases, the total amount of Component 1inside its dissolvable region, and the total mass of the fluid mixture. It is clear that all the quantities shown inFig.13 have little change as time goes on, and we’ve checked that the changes of them with respect to their ini-tial values are below 10 − , which matches our analysis in Section 3.2 and in [21, 23]. The results demonstratethat the proposed scheme conserves the mass of individual pure phases, i.e., (cid:110) ρ φp (cid:80) i,j [ χ p ∆Ω] i,j (cid:111) Np =1 , con-26igure 11: Profiles of the 3 phases in the validation of Galilean invariance. a) Profiles of Phase 1 at x = 0 . x = 0 .
75. c) Profiles of Phase 3 at x = 0 . x = 0 .
5. b) Profiles of Component 2 at x = 0 . (cid:110)(cid:80) i,j [ χ Mp C p ∆Ω] i,j (cid:111) Mp =1 , and, thus, con-serves the mass of the fluid mixture, i.e., (cid:80) i,j [ ρ ∆Ω] i,j = (cid:80) Np =1 ρ φp (cid:80) i,j [ χ p ∆Ω] i,j + (cid:80) Mp =1 ρ Cp (cid:80) i,j [ χ Mp C p ∆Ω] i,j from Eq.(32), in the discrete level.Fig.14 shows the time histories of the changes of the momentum. Using the conservative method, themomentum is conserved at the discrete level. On the other hand, using the balanced-force method doesn’t27igure 14: Time histories of the momentum in the horizontal shear layer. a) Time histories of the x component of the momentum. b) Time histories of the y component of the momentum.Figure 15: Time histories of the energies in the horizontal shear layer. a) Time histories of the kinetic energy,free energy, component energy, and total energy. b) Time histories of the component energy.conserve the momentum exactly. However, the change of the momentum is small, which is less than 0 .
25% ofits initial value in this case. Therefore, the momentum is essentially conserved by the balanced-force method.The obtained results are consistent with our previous work in [24, 22, 21, 23] for multiphase flows withoutcomponents, and they demonstrate that including components doesn’t change the property of momentumconservation of the scheme.Fig.15 shows the time histories of the kinetic energy E K , the free energy E F , the component energy E C and the total energy E T = E K + E F + E C . Both the balanced-force method and the conservative methodgive similar results. In this case, the contribution from the component energy is negligibly small. We canobserve from Fig.15 a) that the kinetic energy is decaying due to both the viscous effect and the energytransfer to the free energy. Therefore, the free energy is increased. The total energy is always decaying,demonstrating that the energy law Eq.(44) is honored at the discrete level by the scheme. Fig.15 b) showsthat the component energy is also decaying all the time. We consider a liquid drop, initially in the air, falling down into another liquid at the bottom, and graduallymixing with the bottom liquid. This is a three-phase case where the liquid drop (Phase 01) is miscible withthe other bottom liquid (Phase 02) while the air (Phase 03) is not miscible with either of the liquids (Phases281 or 02). The densities and viscosities of Phases 01, 02, and 03 are ρ = 3000kg / m , µ = 1 . × − Pa · s, ρ = 1000kg / m , µ = 1 × − Pa · s, ρ = 1kg / m , and µ = 2 × − Pa · s. The surface tension betweenthe immiscible pairs is 0 . / s and the diffusivity between the miscible pair is 1 × − m / s. The gravityis pointing downward with a magnitude 9 . / s .As mentioned in Section 2.6, this three-phase case can be turned into a 2-phase-1-component setup,where Phase 01 is represented by Phase 1 with 1mol / m Component 1, Phase 02 is represented by purePhase 1, Phase 03 is represented by pure Phase 2, and Component 1 is only dissolvable in Phase 1, i.e., { I Mp,q } = { , } . As a result, we have ρ φ = 1000kg / m , µ φ = 1 × − Pa · s, ρ φ = 1kg / m , µ φ = 2 × − Pa · s, σ , = 0 . / s, ρ C = 2000kg / mol, µ C = 5 × − Pa · s · m / mol, D , = 1 × − m / s. The governingequations are non-dimensionalized by a length scale 0 . / m , an acceleration scale1m / s , and a concentration scale 1mol / m . Consequently, Phase 01 is represented by Phase 1 with unityconcentration of Component 1 in the results reported.The domain considered is [1 ×
1] with periodic boundaries along the x axis and with free-slip boundariesalong the y axis. The domain is discretized by [128 × t = 10 − . Initially,the circular drop of Phase 1 is at (0 . , .
75) with a radius 0 .
15, and there is Component 1 with a homogeneousunity concentration inside it. The bottom of the domain below y = 0 . To demonstrate the capability of the proposed model and scheme, we present the final example that includesthree phases and three components, large density and viscosity ratios, multiphase interfacial tensions, andthe effect of contact angles and moving contact lines.The densities and viscosities of the 3 pure phases are ρ φ = 1000kg / m , µ φ = 10 − Pa · s, ρ φ = 500kg / m , µ φ = 10 − Pa · s, ρ φ = 1kg / m , and µ φ = 2 × − Pa · s. The interfacial tensions are σ , = 0 . / m, σ , =0 . / m, and σ , = 0 . / m, and the gravity is g = { , − . } m / s . The contact angle between Phases29igure 16: Results of the miscible falling drop, from left to right and top to bottom, t = 0 .
00, 0 .
20, 0 . .
27, 0 .
30, 0 .
35, 0 .
40, 0 .
50, 0 .
60, 0 .
70, 0 .
80, 0 .
90, 1 .
00, 1 .
20, 1 .
40, 1 .
60, 1 .
80, 2 .
20, 2 .
60, 3 .
00, 3 .
40, 3 .
80, 4 . .
60, 5 .
00. Left of each panel: configuration of the phases, yellow: Phase 1, white: Phase 2. Right of eachpanel: configuration of Component 1, yellow solid lines: interfaces of Phase 1.Figure 17: Time histories of the total volumes of individual phases and the total amount of Component 1 inits dissolvable region in the miscible falling drop. a) Time histories of the total volumes of individual phases.b) Time histories of the total amount of Component 1 in its dissolvable region.1 and 2 at the right boundary is 135 and the other is 90 . The densities and viscosities of the 3 componentsare ρ C = 500kg / mol, µ C = 5 × − Pa · s · m / mol, ρ C = 100kg / mol, µ C = 1 × − Pa · s · m / mol, ρ C = 1kg / mol, µ C = 1 × − Pa · s · m / mol. Component 1 is only dissolvable in Phase 1 with diffusivity1 × − m / s, Component 2 is dissolvable in both Phases 1 and 2 with diffusivity 5 × − m / s and 5 × − m / s, respectively, and Component 3 is dissolvable in both Phases 1 and 3 with diffusivity 5 × − m / sand 2 × − m / s, respectively. The governing equations are non-dimensionalized with a length scale 0 . / m , an acceleration scale 1m / s , and a concentration scale 1mol / m .The domain considered is [1 ×
1] with no-slip boundaries. The domain is discretized by [128 × t = 10 − . Initially, a drop of Phase 1 at (0 . , .
75) with a radius 0 .
15 is above atank of Phase 1 filling 0 (cid:54) y (cid:54) .
3. A drop of Phase 2 is at (0 . , .
6) with a radius 0 .
1. The rest of thedomain is filled by Phase 3. Components 1 and 2 are homogeneously distributed with unity concentration30nside the drops of Phases 1 and 2, respectively. Component 3 is distributed between x = 0 . x = 0 . contact angle between Phases 1 and 2 at the right wall is more clearlyobservable. As time goes on, the phases gradually reach their equilibrium configurations, and Components1 and 3 keep homogenizing inside their corresponding dissolvable regions. Specifically, Component 1 onlyexists inside Phase 1 and none of it is observed in Phases 2 and 3 from our results. On the other hand,Component 2 is allowed in both Phases 1 and 3, and it is not observed in Phase 2 from the results. Again,we plot the time histories of the total volumes of individual phases and the total amounts of each componentin its corresponding dissolvable region in Fig.19, and they are exactly conserved even though the problemconsidered is highly complicated and dynamical.Finally, the effect of the contact angle is illustrated in Fig.20 by comparing the present case to the casewhere all the contact angles are 90 . We observe that the effect of the contact angle doesn’t change muchof the dynamics of the problem while it becomes important for the equilibrium configuration, so we onlyshow the results at t = 5. It is clear that in the case with contact angles all being 90 , i.e., in Fig.20 b), theinterface of Phases 1 and 3 is finally horizontal, and the drop of Phase 2 looks like a section of an ellipsewhich is symmetric with respect to the interface of Phases 1 and 3. On the other hand, in the case with the135 contact angle between Phases 1 and 2 at the right wall, i.e., in Fig.20 a), the interface between Phases1 and 2 is close to an inclined straight line, and the interface between Phases 1 and 3 is not horizontal.The problem considered is again very challenging. It includes multiple materials and their materialproperties are significantly different. The maximum density ratio is 1600, the maximum viscosity ratio is5000, and the diffusivity of the components crosses two orders of magnitude. It also includes gravitation,interfacial tensions, topological change, large deformation of the interfaces, the effects of contact angles, andmoving contact lines, which are the critical factors for multiphase and multicomponent flows. Nevertheless,our model and scheme are capable of overcoming those challenges and producing physically plausible results. In the present work, we developed a consistent and conservative model and its corresponding scheme formultiphase and multicomponent, or N -phase- M -component, incompressible flows with N (cid:62) M (cid:62) t = 0 .
00, 0 .
20, 0 .
25, 0 .
27, 0 .
30, 0 .
35, 0 .
40, 0 .
50, 0 .
60, 0 .
70, 0 .
80, 0 .
90, 1 .
00, 1 .
20, 1 .
40, 1 .
60, 1 .
80, 2 .
20, 2 . .
00, 3 .
40, 3 .
80, 4 .
20, 4 .
60, 5 .
00. Top-left of each panel: configuration of the phases, blue: Phase 1, yellow:Phase 2, white: Phase 3. Top-right of each panel: configuration of Component 1. Bottom-left: configurationof Component 2. Bottom-right: configuration of Component 3. In the top-right, bottom-left, and bottomright of each panel, blue solid lines: interfaces of Phase 1, yellow solid lines: interfaces of Phase 2.phases and components have their own densities and viscosities. Each pair of phases has a surface tensionand each pair of a phase and a component has a diffusivity. At a wall boundary, every two phases have acontact angle.The model is developed based on the multiphase Phase-Field model including the contact angle boundarycondition [11, 12, 21], the diffuse domain approach [36], and the consistency analysis [24] on the consistencyconditions proposed for multiphase and multicomponent flows. The locations of individual phases are rep-resented by a set of order parameters, which is the volume fraction contrasts in the present model, and thesharp interfaces are replaced by interfacial regions inside which there are thermodynamical compression anddiffusion to preserve the thickness of the region. The Phase-Field model also eliminates any generation oflocal void or overfilling, i.e., the summation of the volume fractions from the Phase-Field model is unityeverywhere. Each component has its own concentration governed by the convection-diffusion equation inthe phases that it is dissolvable in, and its flux at the phase boundaries is either zero if it is not dissolvablein the neighboring phases or continuous otherwise. The concentrations represent the local amounts of the32igure 19: Time histories of the total volumes of individual phases and the total amounts of each componentin its corresponding dissolvable region in the falling drops with moving contact lines. a) Time histories ofthe total volumes of individual phases. b) Time histories of the total amounts of each component in itscorresponding dissolvable region.Figure 20: Comparison with different contact angle set-ups at t = 5 in the falling drop with moving contactlines. Left: the contact angle between Phases 1 and 2 is 135 at the right wall. Right: the contact anglebetween Phases 1 and 2 is 90 at the right wall. Top-left of each panel: configuration of the phases,blue: Phase 1, yellow: Phase 2, white: Phase 3. Top-right of each panel: configuration of Component1. Bottom-left: configuration of Component 2. Bottom-right: configuration of Component 3. In the top-right, bottom-left, and bottom-right of each panel, blue solid lines: interfaces of Phase 1, yellow solid lines:interfaces of Phase 2.components and are generally interpreted as the “molar concentrations”, while they can be considered as the“volume fractions” in some special problems. Since all the phases are evolving and deforming, it is challeng-ing to solve the equation and assign the boundary condition of the components. We apply the diffuse domainapproach where the original equation defined in a complex domain along with the boundary condition isreplaced by an equivalent equation defined in the whole domain of interest, with the help of the indicatorfunctions of the phases and phase boundaries. Then, the indicator functions are replaced by the smoothedones, which makes the volume fractions from the Phase-Field model a directly available candidate. The33omponent equation is finalized by incorporating the thermodynamical effect from the Phase-Field model sothat the volume fraction equation derived from the component equation is consistent with the one from thePhase-Field model. We call such a consistency condition, the consistency of volume fraction conservation.This consistency condition is not realized in the original diffuse domain approach [36] and it plays an essentialrole in the energy law and Galilean invariance of the model. The density and viscosity of the fluid mixtureare computed using the volume fractions of the phases and concentrations of the components, and the motionof the fluid mixture is governed by the momentum equation. The consistency of mass conservation and theconsistency of mass and momentum transport are applied to ensure the physical coupling between the massconservation equation and the momentum equation. These consistency conditions are also essential for theenergy law and Galilean invariance of the model. Finally, the analysis of the consistency of reduction isperformed to ensure that the proposed model is not allowed to generate fictitious phases or components. Aphysical energy law is satisfied with the proposed model. We show that, without external input, the totalenergy, which includes free energy, component energy, and kinetic energy, is not increasing as time goeson. Three factors contribute to the decay of the total energy. The first one is due to the thermodynamicalnon-equilibrium in the interfacial regions from the Phase-Field model. The second one is the inhomogeneityof the components in their dissolvable regions from the component equation. The last one is the viscouseffect from the momentum equation. It can also be shown that the proposed model satisfies the Galileaninvariance.Consequently, the proposed model is consistent and conservative. It conserves the mass of individualpure phases, the amount of each component in its dissolvable region, and thus the mass of the fluid mixture,and the momentum of the multiphase and multicomponent flows. It ensures that the summation of thevolume fractions from the Phase-Field model is unity everywhere so that there is no local void or overfilling.It satisfies a physical energy law and it is Galilean invariant. It satisfies all the proposed consistencyconditions, which are the consistency of reduction, the consistency of volume fraction conservation, theconsistency of mass conservation, and the consistency of mass and momentum transport. It should benoted that the consistency conditions play a critical role in avoiding fictitious phases and components, inderiving the energy law, and in proving the Galilean invariance of the model. In addition to multiphaseand multicomponent problems, the proposed model is also applicable to some multiphase problems wherethe miscibilities between each pair of phases are different. The proposed model is also flexible to differentcircumstances of cross-interface transport of a component which is dissolvable in both sides of the interface.Specifically, under different setup, a component, which is dissolvable in two neighboring phases, is either ableto cross the phase interface with continuous flux at the interface or not allowed to cross the phase interface,i.e., with zero flux at the interface. More importantly, the present study proposes a general frameworkthat physically connects the dynamics of the phases, the components, and the fluid flows, using the proposedconsistency conditions. Therefore, it suits to many models for the order parameters or components. Changingthose models is equivalent to updating the definitions of the Phase-Field flux and component flux in thepresent framework and this doesn’t affect the Galilean invariance of and the kinetic energy equation derivedfrom the momentum equation.A few assumptions have been made in the proposed model. First, we have assumed that the concentrationsof the components are governed by the convection-diffusion equation. This is not a strong assumption. Theframework we proposed is flexible to incorporate any other physical model that best describes the dynamicsof the components. Second, we only consider the Neumann-type boundary condition of components atphase interfaces since it suits many applications. The Dirichlet- or Robin-type boundary condition can alsobe incorporated inside the framework of the diffuse domain approach, see [36, 69]. Third, we have notconsidered the effect of components on phase interfaces. For example, the appearance of the components atthe phase interfaces can change the surface tensions as well as the contact angles, and therefore introducethe Marangoni effect. This will be an interesting extension of the present model in future work. Last,the validity of the proposed model relies on the assumption of diluteness, although this assumption canbe removed in a subset of problems where the components dissolved in a specific phase share an identicaldiffusion coefficient in that phase and cross-interface transports of components are not allowed. Completeremoval of this assumption needs to incorporate the effect of components on changing phase volumes andwill be another direction to improve the proposed model.The corresponding numerical scheme is developed for the proposed model. The scheme is formally2nd-order accurate in both time and space. The numerical solution of the model also converges to the34harp-interface one, which is validated by the convergence tests in the present work and in our previouswork [24, 22, 21, 23] for the multiphase incompressible flow model. The scheme preserves the properties ofthe model. It can be shown that the scheme conserves the mass of individual pure phases, the amount ofeach component in its dissolvable region, and therefore the mass of the fluid mixture. The momentum isexactly conserved by the scheme using the conservative method for the interfacial force, while it is essentiallyconserved with the balanced-force method. All the consistency conditions, i.e., the consistency of reduction,the consistency of volume fraction conservation, the consistency of mass conservation, and the consistencyof mass and momentum transport, are shown to be satisfied by the scheme. This is very important for theproblems including large density ratios, see the analysis in [24]. The scheme also ensures that the Phase-Fieldsolutions are always in their physical interval, no fictitious phases or components are numerically generated,and the summation of volume fractions is always unity. The numeral solution also preserves the Galileaninvariance and the energy law, which is demonstrated by numerical experiments. In addition, the capabilityof the model for different circumstances of cross-interface transport of a component which is dissolvable inboth sides of the interface is demonstrated. The effect of the newly proposed consistency of volume fractionconservation on eliminating unphysical fluctuations of the components around phase interfaces is illustrated.Finally, the proposed model and scheme are successfully applied to two challenging problems, includingmultiple materials and various miscibility relations among the materials, and having many critical factorsin multiphase and multicomponent flows, i.e., significant differences of the material properties, gravitation,interfacial tensions, topological change, large deformation of the interfaces, and the effects of contact anglesand moving contact lines. In summary, the proposed model and scheme are general, effective, and robustfor multiphase and multicomponent incompressible flows. 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) [62],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 ScienceFoundation (DMS-1555072 and DMS-1736364).
References [1] D.M. Anderson, G.B. McFadden, and A.A. Wheeler. Diffuse-interface methods in fluid mechanics.
Annu. Rev. Fluid Mech. , 30:139–165, 1998.[2] V.E. Badalassi, H.D. Ceniceros, and S. Banerjee. Computation of multiphase systems with phase fieldmodels.
J. Comput. Phys. , 190:371–397, 2003.[3] F. Boyer and C. Lapuerta. Study of a three component cahn-hilliard flow model.
ESAIM: MathematicalModelling and Numerical Analysis , 40(4):653–687, 2006.[4] F. Boyer, C. Lapuerta, S. Minjeaud, B. Piar, and M. Quintard. Cahn-hilliard/navier-stokes model forthe simulation of three-phase flows.
Transport in Porous Media , 82(3):463–483, 2010.[5] F. Boyer and S. Minjeaud. Hierarchy of consistent n-component cahn–hilliard systems.
Math. ModelsMethods Appl. Sci. , 24:2885–292, 2014.[6] J.U. Brackbill, D.B. Kothe, and C. Zemach. A continuum method for modeling surface tension.
J.Comput. Phys. , 100:335–354, 1992.[7] 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.[8] S. Dong. On imposing dynamic contact-angle boundary conditions for wall-bounded liquid-gas flows.
Comput. Methods Appl. Mech. Engrg. , 247-248:179–200, 2012.359] S. Dong. An efficient algorithm for incompressible n-phase flows.
J. Comput. Phys. , 276:691–728, 2014.[10] S. Dong. Physical formulation and numerical algorithm for simulating n immiscible incompressible fluidsinvolving general order parameters.
J. Comput. Phys. , 836:98–128, 2015.[11] S. Dong. Wall-bounded multiphase flows of nimmiscible incompressible fluids: Consistency and contact-angle boundary condition.
J. Comput. Phys. , 338:21–67, 2017.[12] S. Dong. Multiphase flows of n immiscible incompressible fluids: A reduction-consistent andthermodynamically-consistent formulation and associated algorithm.
J. Comput. Phys. , 361:1–49, 2018.[13] 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.[14] J.H. Ferziger and M. Peric.
Computational Methods for Fluid Dynamics 3rd Edition . Springer, 2002.[15] M.M. Francois. Recent numerical and algorithmic advances within the volume tracking framework formodeling interfacial flows.
Procedia IUTAM , 15:270–277, 2015.[16] 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.[17] F. Gibou, R. Fedkiw, and S. Osher. A review of level-set methods and some recent applications.
J.Comput. Phys. , 353:82–109, 2018.[18] F. Gibou, D. Hyde, and R. Fedkiw. Sharp interface approaches and deep learning techniques formultiphase flows.
J. Comput. Phys. , 380:4420463, 2019.[19] F. Giussani, F. Piscaglia, G. Saez-Mischlich, and J. H`elie. A three-phase vof solver for the simulationof in-nozzle cavitation effects on liquid atomization.
J. Comput. Phys. , 406:109068, 2020.[20] 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.[21] Z. Huang, G. Lin, and A.M. Ardekani. A consistent and conservative phase-field method for multiphaseincompressible flows. arXiv:2010.01099 [physics.comp-ph] , 2020.[22] Z. Huang, G. Lin, and A.M. Ardekani. Consistent and conservative scheme for incompressible two-phaseflows using the conservative allen-cahn model.
J. Comput. Phys. , 420:109718, 2020.[23] 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 [physics.comp-ph] ,2020.[24] Z. Huang, G. Lin, and A.M. Ardekani. Consistent, essentially conservative and balanced-force phase-fieldmethod to model incompressible two-phase flows.
J. Comput. Phys. , 406:109192, 2020.[25] S. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, and L. Tobiska. Quantitativebenchmark computations of two-dimensional bubble dynamics.
Int. J. Numer. Methods. Fluids , 60:1259–1288, 2009.[26] 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.[27] D. Jacqmin. Calculation of two-phase navier-stokes flows using phase-field modeling.
J. Comput. Phys. ,155:96–127, 1999.[28] G-S Jiang and C-W Shu. Efficient implementation of weighted eno schemes.
J. Comput. Phys. ,126:202–228, 1996. 3629] J. Kim. Phase field computations for ternary fluid flows.
Comput. Methods Appl. Mech. Engre. ,196:4779–4788, 2007.[30] J. Kim. A generalized continuous surface tension force formulation for phase-field models for multi-component immiscible fluid flows.
Comput. Methods Appl. Mech. Engre. , 198:3105–3112, 2009.[31] J. Kim. Phase-field models for multi-component fluid flows.
Commun. Comput. Phys. , 12:613–661,2012.[32] J. Kim and H.G. Lee. A new conservative vector-valued allen-cahn equation and its fast numericalmethod.
Comput. Phys. Commun. , 221:102–108, 2017.[33] J. Kim and J. Lowengrub. Phase field modeling and simulation of three-phase flows.
Interfaces FreeBound. , 7:435–466, 2005.[34] 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.[35] H.G. Lee and J. Kim. An efficient numerical method for simulating multiphase flows using a diffuseinterface model.
Physica A , 423:33–50, 2015.[36] X. Li, J. Lowengrub, A. Ratz, and A. Voigt. Solving pdes in complex geometries: A diffuse domainapproach.
Commun. Math. Sci. , 1:81–107, 2009.[37] F. Losasso, T. Shinar, A. Selle, and R. Fedkiw. Multiple interacting liquids.
ACM Transactions onGraphics (TOG) , 25(3):812–819, 2006.[38] E. Olsson and G. Kreiss. A conservative level set method for two phase flow.
J. Comput. Phys. ,210:225–246, 2005.[39] E. Olsson, G. Kreiss, and S. Zahedi. A conservative level set method for two phase flow ii.
J. Comput.Phys. , 225:785–807, 2007.[40] S. Osher and A.J. Sethian. Fronts propagating with curvature-dependent speed: Algorithms based onhamilton-jacobi formulations.
J. Comput. Phys. , 79:12–49, 1988.[41] M. Owkes and O. Desjardins. A mass and momentum conserving unsplit semi-lagrangian frameworkfor simulating multiphase flows.
J. Comput. Phys. , 332:21–46, 2017.[42] A. Pathak and M. Raessi. A three-dimensional volume-of-fluid method for reconstructing and advectingthree-material interfaces forming contact lines.
J. Comput. Phys. , 307:550–573, 2016.[43] S. Popinet. An accurate adaptive solver for surface-tension-driven interfacial flows.
J. Comput. Phys. ,228:5838–5866, 2009.[44] S. Popinet. Numerical models for surface tension.
Annu. Rev. Fluid Mech. , 50:49–75, 2018.[45] A. Prosperetti and G. Tryggvason.
Computational Methods for Multiphase Flow . Cambridge UniversityPress, 2007.[46] 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.[47] Abadi R.H.H., M.H. Rahimian, and A. Fakhari. Conservative phase-field lattice-boltzmann model forternary fluids.
J. Comput. Phys. , 374:668–691, 2018.[48] N. Scapin, P. Costa, and L. Brandt. A volume-of-fluid method for interface-resolved simulations ofphase-changing two-fluid flows.
J. Comput. Phys. , 407:109251, 2020.[49] R. Scardovelli and S. Zaleski. Direct numerical simulation of free-surface and interfacial flow.
Annu.Rev. Fluid Mech. , 31:567–603, 1999. 3750] S.P. Schofield, R.V. Garimella, M.M. Francois, and R. Loubere. A second-order accurate material-order-independent interface reconstruction technique for multi-material flow simulations.
J. Comput. Phys. ,228:731–745, 2009.[51] S.P. Schofield, Christon M.A., V. Dyadechko, R.V. Garimella, R.B. Lowrie, and B.K. Swartz. Multi-material incompressible flow simulation using the moment-of-fluid method.
Int. J. Numer. Meth. Fluids ,63:931–952, 2010.[52] J.A. Sethian and P. Smereka. Level set method for fluid interfaces.
Annu. Rev. Fluid Mech. , 35:341–372,2003.[53] 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.[54] Y. Shi, G.H. Tang, L.H. Cheng, and H.Q. Shuang. An improve d phase-field-base d lattice boltzmannmodel for droplet dynamics with soluble surfactant.
Comput. Fluids , 179:508–520, 2019.[55] K.A. Smith, F.J. Solis, and D.L. Chopp. A projection method for motion of triple junctions by levelsets.
Interfaces Free Bound. , 4:263–276, 2002.[56] G. Soligo, A. Roccon, and A. Soldati. Coalescence of surfactant-laden drops by phase field method.
J.Comput. Phys. , 376:1292–1311, 2019.[57] D.P. Starinshak, S. Karni, and P.L. Roe. A new level set model for multimaterial flows.
Interface andfree boundaries , 4:263–276, 2002.[58] 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.[59] M. Sussman, K.M. Smith, M.Y. Hussaini, M. Ohta, and R. Zhi-Wei. A sharp interface method forincompressible two-phase flows.
Journal of computational physics , 221(2):469–505, 2007.[60] ML Szulczewski and R Juanes. The evolution of miscible gravity currents in horizontal porous layers.
Journal of Fluid Mechanics , 719:82–96, 2013.[61] K.E. Teigen, P. Song, J. Lowengrub, and A. Voigt. A diffuse-interface method for two-phase flows withsoluble surfactants.
J. Comput. Phys. , 230:375–393, 2011.[62] 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.[63] 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.[64] G. Tryggvason, R. Scardovelli, and S. Zaleski.
Direct Numerical Simulations of Gas-Liquid MultiphaseFlows . Cambridge University Press, 2011.[65] S.O. Unverdi and G. Tryggvason. A front-tracking method for viscous, incompressible, multi-fluid flows.
J. Comput. Phys. , 100:25–37, 1992.[66] S. Wu and J. Xu. Multiphase allen-cahn and cahn-hilliard models and their discretizations with theeffect of pairwise surface tensions.
J. Comput. Phys. , 343:10–32, 2017.[67] 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.[68] 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. 3869] F. Yu, Z. Guo, and J. Lowengrub. Higher-order accurate diffuse-domain methods for partial differen-tial equations with dirichlet boundary conditions in complex, evolving geometries.
J. Comput. Phys. ,406:109174, 2020.[70] C.Y. Zhang, H. Ding, P. Gao, and Y.L. Wu. Diffuse interface simulation of ternary fluids in contactwith solid.
J. Comput. Phys. , 309:37–51, 2016.[71] Q. Zhang and X.P. Wang. Phase field modeling and simulation of three-phase flow on solid surfaces.
J.Comput. Phys. , 319:79–107, 2016.[72] G. Zhu, J. Kou, J. Yao, A. Li, and S. Sun. A phase-field moving contact line model with solublesurfactants.