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

The problem considered in the present work includes multiple phases and components, and we use $N$ ($N\geqslant 1$) to denote the number of phases and $M$ ($M\geqslant 0$) for the number of components.

The $N$ phases are specified by defining a set of order parameters $\{\phi_{p}\}_{p=1}^{N}\in[-1,1]$, which are the volume fraction contrasts of individual phases. At the location of $\phi_{p}=1$, there is only Phase $p$, while $\phi_{p}=-1$ represents the absence of Phase $p$. The volume fractions of individual phases are related to $\{\phi_{p}\}_{p=1}^{N}$ by

$$\chi_{p}=\frac{1+\phi_{p}}{2},1\leqslant p\leqslant N,$$

(1)

so that $\{\chi_{p}\}_{p=1}^{N}\in[0,1]$. Every location should be at least occupied by one of the phases while at most occupied by all of them. As a result, the volume fractions are not independent of each other and should satisfy the summation constraint

$$\sum_{p=1}^{N}\chi_{p}=1,$$

(2)

or equivalently

$$\sum_{p=1}^{N}\phi_{p}=2-N,$$

(3)

everywhere. Each phase has a background fluid called the pure phase, and $\rho_{p}^{\phi}$ and $\mu_{p}^{\phi}$ denote the density and viscosity, respectively, of pure Phase $p$. Every pair of phases is immiscible, and therefore there is a pairwise surface tension between them, denoted by $\sigma_{p,q}$ for Phases $p$ and $q$. The pairwise surface tensions build a $N\times N$ matrix which is symmetric and has zero diagonal. At a wall boundary, there is a pairwise contact angle $\theta_{p,q}^{W}$ of Phases $p$ and $q$, measured on the side of Phase $p$. Again, the pairwise contact angles build a $N\times N$ matrix whose entries satisfy $\theta_{p,q}^{W}+\theta_{q,p}^{W}=\pi$. As already mentioned, the Marangoni effect of the components is assumed to be negligible, and consequently, the pairwise surface tensions and contact angles do not change due to the appearance of the components.

The $M$ components are represented by a set of concentrations $\{C_{p}\}_{p=1}^{M}$. Individual components have their densities $\{\rho_{p}^{C}\}_{p=1}^{M}$ and viscosities $\{\mu_{p}^{C}\}_{p=1}^{M}$. Each component can be dissolved inside different phases. On the other hand, there can be multiple components inside a specific phase. To indicate the dissolvabilities among the components and the phases, the dissolvability matrix with dimension $M\times N$ is defined as

$$I^{M}_{p,q}=\left\{\begin{array}[]{l}1,\textrm{if Component $p$ is dissolvable in Phase $q$},\\ 0,\textrm{else}.\end{array}\right.$$

(4)

The diffusion coefficient of Component $p$ in Phase $q$ is denoted by $D_{p,q}$. Noted that values of $\{C_{p}\}_{p=1}^{M}$ are meaningful only inside the corresponding phases they dissolve in. To obtain a meaningful value in the entire domain of interest, one needs to use, e.g., $I^{M}_{p,q}\chi_{q}C_{p}$, to represent the amount of Component $p$ in Phase $q$. As pointed out in SzulczewskiJuanes2013 , there are two commonly used ways to interpret concentration. The first case is to quantify the amount of a dissolved solute of a solution such as the “molar concentration” of salt in a salt water solution. The second case is to quantify the composition of miscible fluids such as the “volume fraction” of ethanol 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 “molar concentration”. We further assume that inside each phase the “solution” compound of the pure phase (as the “solvent”) and the components (as the “solutes”) dissolved in that phase is dilute. In other words, in a specific phase, the volume fractions of the components dissolved in it are negligibly small compared to the volume fraction of the pure phase or the background fluid of that phase. Therefore, the volume of each phase will not be changed by either the appearance or cross-interface transport of the components. This assumption of diluteness is also consistent with neglecting the Marangoni effect of the components. It should be 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 all the components dissolved in it and the cross-interface transport of those components is prohibited. Although the assumption of diluteness can be removed in these problems, one has to be cautious about the assumption of neglecting the Marangoni effect of the components when applying the proposed model to these problems. In the present work, we loosely call $\{\rho_{p}^{C}\}_{p=1}^{M}$ and $\{\mu_{p}^{C}\}_{p=1}^{M}$ densities and viscosities, respectively, regardless of which interpretation of the concentrations is employed. The actual meaning of $\{\rho_{p}^{C}\}_{p=1}^{M}$ depends on the concrete definition of the concentrations. For example, $\rho_{p}^{C}$ is the molar mass of Component $p$ when the concentrations are “molar concentration”, while it is the density of pure Component $p$ minus the density of the pure phase the component dissolves in when the concentrations mean “volume fraction”. The same works for $\{\mu_{p}^{C}\}_{p=1}^{M}$.

The densities of the pure phases and components are all constant, and phase changes are not considered in the present work, in addition to the aforementioned assumptions. As a result, the velocity of the flow $\mathbf{u}$ is divergence-free, i.e.,

$$\nabla\cdot\mathbf{u}=0.$$

(5)

The governing equations of the $N$-phase-$M$-component model are described in detail, which include the multiphase Phase-Field model, the diffuse domain approach for the components, the material properties, and the momentum equation.

The consistency analysis is of great importance for building up a physical multiphase model, as shown in our previous studies of Phase-Field modeling for two- and multi-phase flows (Huangetal2020, ; Huangetal2020N, ; Huangetal2020CAC, ). Specifically, the consistency 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^{\prime}$-phase-$M^{\prime}$-component system ($1\leqslant N^{\prime}\leqslant N-1$, $0\leqslant M^{\prime}\leqslant M-1$), when ($N-N^{\prime}$) phases and $(M-M^{\prime})$ components are absent.

• •

  Consistency of volume fraction conservation: The conservation equation of the volume fractions derived from the component equation should be consistent with the one derived from the Phase-Field equation when the component is homogeneous.

• •

  Consistency of mass conservation: The mass conservation equation should be consistent with the one derived from the Phase-Field equation, the component equation, and the density of the fluid mixture. The mass flux $\mathbf{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 should be computed as a tensor product between the mass flux and the velocity, where the mass flux should be identical to the one in the mass conservation equation.

The consistency of reduction is enriched from the $N$-phase one by adding the component part. It is important because it assures that no fictitious phase or component is generated by the model. Any phases or components that are absent at $t=0$ will not appear at $\forall 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 $\{\phi_{p}\}_{p=1}^{N}$, has been shown in (Dong2018, ; Huangetal2020N, ). The momentum equation is reduction consistent as long as the density $\rho$, the viscosity $\mu$, and the mass flux $\mathbf{m}$ are consistent as well. Thus, our attention is paid on the part related to the components. For a specific $p$ and from Eq.(25), we can obtain $C_{p}\equiv 0,\forall 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 the mass 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^{\prime})$ times so that the $M$-component system reduces to the corresponding $M^{\prime}$-component system ($0\leqslant M^{\prime}\leqslant M-1$). Combining with the analysis in (Huangetal2020N, ), the proposed $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 diffuse domain approach, and has been considered in Section 2.2.2 when building up the component equation Eq.(25). In the original diffuse domain approach, it presumes that the smoothed indicator function has the same sharp-interface dynamics as the exact one Eq.(22), i.e., the phase interface is advected by the flow velocity. 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 but finite. 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 original diffuse domain approach. The consistency of volume fraction conservation aims to remove the discrepancy by incorporating the thermodynamical effects in the Phase-Field model to the component equation. This consistency condition is essential for the model satisfying the energy law and the Galilean invariance, see Sections 2.4 and 2.5.

The consistency of mass conservation and the consistency of mass and momentum transport are identical to those for two- and multi-phase flows (Huangetal2020, ; Huangetal2020N, ; Huangetal2020CAC, ). Based on the definitions of the density of the fluid mixture Eq.(32) and the consistent mass flux Eq.(35), it can be shown that the density is governed by

$$\frac{\partial\rho}{\partial t}+\nabla\cdot\mathbf{m}=\sum_{p=1}^{N}\rho_{p}^{\phi}\underbrace{\left(\frac{\partial\chi_{p}}{\partial t}+\nabla\cdot\mathbf{m}_{\chi_{p}}\right)}_{0}+\sum_{p=1}^{M}\rho_{p}^{C}\underbrace{\left(\frac{\partial(\chi_{p}^{M}C_{p})}{\partial t}+\nabla\cdot\mathbf{m}_{C_{p}}\right)}_{0}=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 the consistent mass flux Eq.(35) has been applied to the inertial term of the momentum equation Eq.(34), the consistency of mass and momentum transport is satisfied as well. These two consistency conditions are critical for problems with large density ratios, as indicated in our previous studies (Huangetal2020, ; Huangetal2020N, ; Huangetal2020CAC, ), since they honor the physical coupling between the mass and momentum. Moreover, they are essential for the model satisfying 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 free energy density, component energy density, and kinetic energy density to be

$$e_{F}=\sum_{p,q=1}^{N}\frac{\lambda_{p,q}}{2}\left(\frac{1}{\eta^{2}}\left(g_{1}(\phi_{p})+g_{1}(\phi_{q})-g_{2}(\phi_{p}+\phi_{q})\right)-\nabla\phi_{p}\cdot\nabla\phi_{q}\right),$$

(38)

$$e_{C}=\sum_{p=1}^{M}\frac{1}{2}\gamma_{C_{p}}\chi_{p}^{M}C_{p}^{2},$$

(39)

$$e_{K}=\frac{1}{2}\rho\mathbf{u}\cdot\mathbf{u},$$

(40)

respectively, where $\{\gamma_{C_{p}}\}_{p=1}^{M}$ are positive dimensional constants so that $e_{C}$ has a unit $[\mathrm{J/m^{3}}]$ and we choose them to be unity without loss of generality. The corresponding energies are the integrals of the energy densities over the domain. It should be noted that the chemical potential of Phase $p$, i.e., Eq.(8), is the functional derivative of the free energy with respect to $\phi_{p}$, i.e., $\xi_{p}=\frac{\delta E_{F}}{\delta\phi_{p}}=\frac{\delta}{\delta\phi_{p}}\int_{\Omega}e_{F}d\Omega$.

After multiplying the Phase-Field equation Eq.(6) with the chemical potential Eq.(8) and then summing over all the phases $p$, we obtain the governing equation for the free energy density, i.e.,

	$\displaystyle\frac{1}{2}\frac{\partial e_{F}}{\partial t}+\sum_{p.q=1}^{N}\frac{\lambda_{p,q}}{4}\nabla\cdot\left(\frac{\partial\phi_{p}}{\partial t}\nabla\phi_{q}+\frac{\partial\phi_{q}}{\partial t}\nabla\phi_{p}\right)+\mathbf{u}\cdot\frac{1}{2}\sum_{p=1}^{N}\xi_{p}\nabla\phi_{p}$		(41)
	$\displaystyle=\frac{1}{2}\sum_{p,q=1}^{N}\nabla\cdot\left(\xi_{p}M_{p,q}^{\phi}\nabla\xi_{q}\right)-\frac{1}{2}\sum_{p,q=1}^{N}M_{p,q}^{\phi}\nabla\xi_{p}\cdot\nabla\xi_{q}.$

After multiplying the component equation Eq.(25) with $\gamma_{C_{p}}C_{p}$ and then summing over all the components $p$, we obtain the governing equation for the component energy density, i.e.,

$$\frac{\partial e_{C}}{\partial t}+\sum_{p=1}^{M}\nabla\cdot\left(\mathbf{m}_{\chi_{p}^{M}}\frac{1}{2}\gamma_{C_{p}}C_{p}^{2}\right)=\sum_{p=1}^{M}\nabla\cdot(\gamma_{C_{p}}D_{p}C_{p}\nabla C_{p})-\sum_{p=1}^{M}\gamma_{C_{p}}D_{p}\nabla C_{p}\cdot\nabla C_{p}.$$

(42)

Eq.(23) has been used in the derivation.

After performing the dot product between $\mathbf{u}$ and the momentum equation Eq.(34), we obtain the governing equation for the kinetic energy density, i.e.,

$\displaystyle\frac{\partial e_{K}}{\partial t}+\nabla\cdot\left(\mathbf{m}\frac{\mathbf{u}\cdot\mathbf{u}}{2}\right)=-\nabla\cdot(\mathbf{u}P)+\nabla\cdot\left(\mu(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})\cdot\mathbf{u}\right)-\frac{1}{2}\mu(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}):(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})+\mathbf{u}\cdot\mathbf{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 the domain boundary, we obtain the energy law for the $N$-phase-$M$-component system, i.e.,

	$\displaystyle\frac{dE_{T}}{dt}=\frac{d}{dt}\left(\frac{1}{2}E_{F}+E_{C}+E_{K}\right)=\frac{d}{dt}\int_{\Omega}\left(\frac{1}{2}e_{F}+e_{C}+e_{K}\right)d\Omega$		(44)
	$\displaystyle=-\int_{\Omega}\frac{1}{2}\sum_{p,q=1}^{N}M_{p,q}^{\phi}\nabla\xi_{p}\cdot\nabla\xi_{q}d\Omega-\int_{\Omega}\sum_{p=1}^{M}\gamma_{C_{p}}D_{p}\nabla C_{p}\cdot\nabla C_{p}d\Omega-\int_{\Omega}\frac{1}{2}\mu(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}):(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})d\Omega.$

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. The right-hand side (RHS) of Eq.(44) shows three factors to reduce the total energy of the multiphase and multicomponent system. The first one comes from the thermodynamical non-equilibrium, contributed by the Phase-Field model. It should be noted that it is effective only in the interfacial region since the mobility Eq.(7) is zero inside the bulk-phase region. The second one results from the inhomogeneity of the components in 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 due to the viscosity of the fluids, contributed from the momentum equation. To include the effect of the contact angles, i.e., the contact angles are other than $\frac{\pi}{2}$, a wall energy functional has to be introduced. However, so far, it is still an open question whether the contact angle boundary condition Eq.(13) from (Dong2017, ) corresponds to a multiphase wall functional.

At the end of this section, we need to emphasize the significance of the consistency conditions in deriving the energy law. The left-hand side (LHS) of Eq.(42) is derived from

	$\displaystyle C_{p}\left(\frac{\partial(\chi_{p}^{M}C_{p})}{\partial t}+\nabla\cdot(\mathbf{m}_{\chi_{p}^{M}}C_{p})\right)=C_{p}\left(\chi_{p}^{M}\frac{\partial C_{p}}{\partial t}+\mathbf{m}_{\chi_{p}^{M}}\cdot\nabla C_{p}\right)+C_{p}^{2}\underbrace{\left(\frac{\partial\chi_{p}^{M}}{\partial t}+\nabla\cdot\mathbf{m}_{\chi_{p}^{M}}\right)}_{0}$		(45)
	$\displaystyle=\chi_{p}^{M}\frac{\partial\frac{1}{2}C_{p}^{2}}{\partial t}+\mathbf{m}_{\chi_{p}^{M}}\cdot\nabla\frac{1}{2}C_{p}^{2}+\frac{C_{p}^{2}}{2}\underbrace{\left(\frac{\partial\chi_{p}^{M}}{\partial t}+\nabla\cdot\mathbf{m}_{\chi_{p}^{M}}\right)}_{0}$
	$\displaystyle=\frac{\partial(\frac{1}{2}\chi_{p}^{M}C_{p}^{2})}{\partial t}+\nabla\cdot\left(\mathbf{m}_{\chi_{p}^{M}}\frac{1}{2}C_{p}^{2}\right).$

Thanks to the consistency of volume fraction conservation, the volume fraction equation, i.e., the terms with under-brace in Eq.(45), is recovered. Similarly, the derivation of the left-hand side (LHS) of Eq.(43) is

$\displaystyle\mathbf{u}\cdot\left(\frac{\partial(\rho\mathbf{u})}{\partial t}+\nabla\cdot(\mathbf{m}\otimes\mathbf{u})\right)=\frac{\partial e_{K}}{\partial t}+\nabla\cdot\left(\mathbf{m}\frac{1}{2}\mathbf{u}\cdot\mathbf{u}\right)+\frac{1}{2}\mathbf{u}\cdot\mathbf{u}\underbrace{\left(\frac{\partial\rho}{\partial t}+\nabla\cdot\mathbf{m}\right)}_{0}=\frac{\partial e_{K}}{\partial t}+\nabla\cdot\left(\mathbf{m}\frac{1}{2}\mathbf{u}\cdot\mathbf{u}\right).$

(46)

The mass conservation equation, i.e., Eq.(37), is recovered in Eq.(46), thanks to the consistency of mass conservation and the consistency of mass and momentum transport. If the consistency conditions are violated, the terms with under-brace in Eq.(45) and Eq.(46) are not zero any more. As a result, some additional terms will 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 homogeneous in their dissolvable region, and all the fluids are inviscid, the total energy of the multiphase and multicomponent system could still change due to the additional terms introduced by the inconsistency. Such a behavior of the total energy is not physically plausible. Thus, the consistency conditions are playing an essential role in the physical behavior of the system.

The proposed $N$-phase-$M$-component model satisfies the Galilean invariance, thanks to the consistency conditions. In this section, we call the reference frame Frame $(\mathbf{x},t)$, and another frame that is moving with respect to the reference frame with a constant velocity $\mathbf{u}_{0}$ Frame $(\mathbf{x}^{\prime},t^{\prime})$. All the variables with ^′ are measured in Frame $(\mathbf{x}^{\prime},t^{\prime})$. The Galilean transform is

$$\mathbf{x}^{\prime}=\mathbf{x}-\mathbf{u}_{0}t,\quad t^{\prime}=t,\quad\mathbf{u^{\prime}}=\mathbf{u}-\mathbf{u}_{0},\quad f^{\prime}=f,\quad\nabla^{\prime}f^{\prime}=\nabla f,\quad\nabla^{\prime 2}f^{\prime}=\nabla^{2}f,\quad\frac{\partial f^{\prime}}{\partial t^{\prime}}=\frac{\partial f}{\partial t}+\mathbf{u}_{0}\cdot\nabla 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 $(\mathbf{x}^{\prime},t^{\prime})$ are

$$\nabla^{\prime}\cdot\mathbf{u}^{\prime}=\nabla\cdot(\mathbf{u}-\mathbf{u}_{0})=0,$$

(48)

$\displaystyle\frac{\partial\phi^{\prime}_{p}}{\partial t^{\prime}}+\nabla^{\prime}\cdot(\mathbf{u}^{\prime}\phi^{\prime}_{p})-\nabla^{\prime}\cdot\left(\sum_{q=1}^{N}M^{\prime\phi}_{p,q}\nabla^{\prime}\xi^{\prime}_{q}\right)=\frac{\partial\phi_{p}}{\partial t}+\mathbf{u}_{0}\cdot\nabla\phi_{p}+\nabla\cdot\left((\mathbf{u}-\mathbf{u}_{0})\phi_{p}\right)-\nabla\cdot\left(\sum_{q=1}^{N}M_{p,q}^{\phi}\nabla\xi_{q}\right)=0,$

(49)

$\displaystyle\mathbf{m}_{\phi^{\prime}_{p}}=\mathbf{u}^{\prime}\phi^{\prime}_{p}-\sum_{q=1}^{N}M^{\prime}_{p,q}\nabla^{\prime}\xi^{\prime}_{q}=(\mathbf{u}-\mathbf{u}_{0})\phi_{p}-\sum_{q=1}^{N}M_{p,q}\nabla\xi_{q}=\mathbf{m}_{\phi_{p}}-\mathbf{u}_{0}\phi_{p},$

(50)

$\displaystyle\mathbf{m}_{\chi^{\prime}_{p}}=\frac{1}{2}(\mathbf{u}^{\prime}+\mathbf{m}_{\phi^{\prime}_{p}})=\frac{1}{2}(\mathbf{u}-\mathbf{u}_{0}+\mathbf{m}_{\phi_{p}}-\mathbf{u}_{0}\phi_{p})=\mathbf{m}_{\chi_{p}}-\mathbf{u}_{0}\chi_{p},$

(51)

$\displaystyle\mathbf{m}_{\chi^{\prime M}_{p}}=\sum_{q=1}^{N}I_{p,q}^{M}\mathbf{m}_{\chi^{\prime}_{q}}=\sum_{q=1}^{M}I_{p,q}^{M}\mathbf{m}_{\chi_{q}}-\mathbf{u}_{0}\sum_{q=1}^{M}I_{p,q}^{M}\chi_{q}=\mathbf{m}_{\chi_{p}^{M}}-\mathbf{u}_{0}\chi_{p}^{M},$

(52)

	$\displaystyle\frac{\partial(\chi^{\prime M}_{p}C^{\prime}_{p})}{\partial t^{\prime}}+\nabla^{\prime}\cdot(\mathbf{m}_{\chi^{\prime M}_{p}}C^{\prime}_{p})-\nabla^{\prime}\cdot(D^{\prime}_{p}\nabla^{\prime}C^{\prime}_{p})=$		(53)
	$\displaystyle\frac{\partial(\chi_{p}^{M}C_{p})}{\partial t}+\mathbf{u}_{0}\cdot\nabla(\chi_{p}^{M}C_{p})+\nabla\cdot\left((\mathbf{m}_{\chi_{p}^{M}}-\mathbf{u}_{0}\chi_{p}^{M})C_{p}\right)-\nabla\cdot(D_{p}\nabla C_{p})=0,$

$\displaystyle\mathbf{m}_{C^{\prime}_{p}}=\mathbf{m}_{\chi^{\prime M}_{p}}C^{\prime}_{p}-D^{\prime}_{p}\nabla^{\prime}C^{\prime}_{p}=(\mathbf{m}_{\chi_{p}^{M}}-\mathbf{u}_{0}\chi_{p}^{M})C_{p}-D_{p}\nabla C_{p}=\mathbf{m}_{C_{p}}-\mathbf{u}_{0}\chi_{p}^{M}C_{p},$

(54)

$\displaystyle\mathbf{m}^{\prime}=\sum_{p=1}^{N}\rho_{p}^{\phi}\mathbf{m}_{\chi^{\prime}_{p}}+\sum_{p=1}^{M}\rho_{p}^{C}\mathbf{m}_{C^{\prime}_{p}}=\sum_{p=1}^{N}\rho_{p}^{\phi}(\mathbf{m}_{\chi_{p}}-\mathbf{u}_{0}\chi_{p})+\sum_{p=1}^{M}\rho_{p}^{C}(\mathbf{m}_{C_{p}}-\mathbf{u}_{0}\chi_{p}^{M}C_{p})=\mathbf{m}-\mathbf{u}_{0}\rho,$

(55)

$\displaystyle\frac{\partial\rho^{\prime}}{\partial t^{\prime}}+\nabla^{\prime}\cdot\mathbf{m}^{\prime}=\frac{\partial\rho}{\partial t}+\mathbf{u}_{0}\cdot\nabla\rho+\nabla\cdot(\mathbf{m}-\mathbf{u}_{0}\rho)=0,$

(56)

	$\displaystyle\frac{\partial(\rho^{\prime}\mathbf{u}^{\prime})}{\partial t^{\prime}}+\nabla^{\prime}\cdot(\mathbf{m}^{\prime}\otimes\mathbf{u}^{\prime})+\nabla^{\prime}P^{\prime}-\nabla^{\prime}\cdot[\mu^{\prime}(\nabla^{\prime}\mathbf{u}^{\prime}+\nabla^{\prime}\mathbf{u}^{\prime T})]-\rho^{\prime}\mathbf{g}-\mathbf{f}^{\prime}_{s}=$		(57)
	$\displaystyle\frac{\partial(\rho(\mathbf{u}-\mathbf{u}_{0}))}{\partial t}+\mathbf{u}_{0}\cdot\nabla(\rho(\mathbf{u}-\mathbf{u}_{0}))+\nabla\cdot\left((\mathbf{m}-\rho\mathbf{u}_{0})\otimes(\mathbf{u}-\mathbf{u}_{0})\right)$
	$\displaystyle+\nabla P-\nabla\cdot[\mu(\nabla(\mathbf{u}-\mathbf{u}_{0})+\nabla(\mathbf{u}-\mathbf{u}_{0})^{T})]-\rho\mathbf{g}-\mathbf{f}_{s}=$
	$\displaystyle\frac{\partial(\rho\mathbf{u})}{\partial t}+\nabla\cdot(\mathbf{m}\otimes\mathbf{u})-\mathbf{u}_{0}\underbrace{\left(\frac{\partial\rho}{\partial t}+\nabla\cdot\mathbf{m}\right)}_{0}+\nabla P-\nabla\cdot[\mu(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})]-\rho\mathbf{g}-\mathbf{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 $(\mathbf{x}^{\prime},t^{\prime})$ are identical to their correspondences, i.e., Eq.(6), Eq.(25), Eq.(37), Eq.(34) and Eq.(48), in Frame $(\mathbf{x},t)$. Thus the Galilean invariance is satisfied. It should be noted that the consistency of mass conservation and the consistency of mass and momentum transport play critical roles in the Galilean invariance of the momentum equation, see the under-braced term in Eq.(57). If either of the consistency conditions is violated, the under-braced term in 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 momentum equation. Consider the case where $\mathbf{u}^{\prime}=\mathbf{0}$, the divergence-free condition is satisfied and the momentum equation in Frame $(\mathbf{x}^{\prime},t^{\prime})$ requires that

$$-\nabla^{\prime}P^{\prime}+\rho^{\prime}\mathbf{g}+\mathbf{f}^{\prime}_{s}=\mathbf{0},$$

(58)

which corresponds to the mechanical equilibrium. From the Galilean transform, we have the same mechanical equilibrium, i.e.,

$$-\nabla P+\rho\mathbf{g}+\mathbf{f}_{s}=\mathbf{0},$$

(59)

and $\mathbf{u}=\mathbf{u}_{0}$ in Frame $(\mathbf{x},t)$. So the divergence-free condition is true in Frame $(\mathbf{x},t)$. The momentum equation in Frame $(\mathbf{x},t)$ becomes

$$\frac{\partial(\rho\mathbf{u}_{0})}{\partial t}+\nabla\cdot(\mathbf{m}\otimes\mathbf{u}_{0})=\nabla\cdot[\mu(\nabla\mathbf{u}_{0}+\nabla\mathbf{u}_{0}^{T})]$$

(60)

The right-hand side (RHS) is zero due to $\nabla\mathbf{u}_{0}=\mathbf{0}$. The left-hand side(LHS) is again zero since the consistency of mass conservation and the consistency of mass and momentum transport recover the mass conservation equation Eq.(37). Consequently, the homogeneous velocity $\mathbf{u}_{0}$ is an admissible solution.

A consistent and conservative $N$-phase-$M$-component flow model is proposed, which allows different densities and 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 consistency conditions, conserves the mass of individual pure phases, conserves the amount of individual components inside 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 dependent on 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 of phases are different. For example, to model a three-phase system where Phases 01 and 02 are miscible while Phases 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., $1\mathrm{mol/m^{3}}$, of Component 1 represents Phase 02, and Phase 2 represents Phase 03. We can set $I_{1,1}^{M}=1$ and $I_{1,2}^{M}=0$. As a result, Phase 1, with/without Component 1, is immiscible with Phase 2, and the immisciblities between Phases 01 and 03 and between Phases 02 and 03 are properly modeled. Component 1 is only dissolvable in 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 models the miscible behavior of Phases 01 and 02. However, the proposed model is not ready for the problems where the miscible phases have significantly different surface tensions with respect to a phase that they are immiscible 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 volumes due to the transport of components has been considered.