Advances of Phase-Field Model in the Numerical Simulation of Multiphase Flows: A Review

Abstract

The phase-field model (PFM) is gaining increasing attention in the application of multiphase flows due to its advantages, in which the phase interface is treated as a narrow layer and phase parameters change smoothly and continually at this thin layer. Thus, the construction or tracking of the phase interface can be avoided, and the bulk phase and phase interface can be simulated integrally. PFM provides a useful alternative that does not suffer from problems with either the mass conservation or the accurate computation of surface tension. In this paper, the state of the art of PFM in the numerical modeling and simulation of multiphase flows is comprehensively reviewed. Starting with a brief description of historical developments in the PFM, we continue to take a tour into the basic concepts, fundamental theory, and mathematical models. Then, the commonly used numerical schemes and algorithms for solving the governing systems of PFM in the application of multiphase flows are presented. The various applications and representative results, especially in non-match density scenarios of multiphase flows, are reviewed. The primary challenges and research focus of PFM are analyzed and summarized as well. This review is expected to provide a valuable reference for PFM in the application of multiphase flows.

1. Introduction

The difference and diversity of material properties contribute to the essence of the world. Phase, as one of the material properties, defines the part with uniform physical and thermal properties in a concerned thermodynamic system that can be separated naturally by an interface [1]. Multiphase flow has found a wide range of applications in various aspects of engineering and science fields, such as petroleum engineering, mechanical engineering, and chemical engineering, etc. [1,2,3]. However, whether multiphase flow is favorable or not depends on a deep understanding of its physical essence. Therefore, the study of multiphase flow has enormous scientific significance and engineering application value.

In recent years, the development of production technology and the complexity of engineering problems have motivated diverse researches in multiphase flows. The commonly used approaches for multiphase flow investigation include theoretical analysis, experimental study, and numerical simulation [2], etc., among which the theoretical analysis is capable of accurately providing the mathematical relationship between various phase parameters. However, it is only applicable to simple multiphase flow problems. The experimental study is close to engineering practice, but it is usually limited by the high cost, long experimental period, and huge difficulties to be conducted in some special scenarios. In addition, the semi-empirical or empirical formulas obtained from experimental research are not applicable to various applications of multiphase flow. In addition, the influencing factors are usually difficult to separately investigate in the experimental study. Comparatively, the numerical simulation which is performed based on mathematical models can be used as an essential supplement for the theoretical analysis and experimental study. Moreover, the cost of numerical simulation is relatively low and the research period is relatively short. The obtained simulation data can help reveal the underlying mechanism behind the ‘flow pattern’ in detail. Therefore, the numerical simulation has been a highly appreciated approach in the study of multiphase flows.

Significant advances have been witnessed in the numerical modeling and simulation techniques for multiphase flows in past decades. Fuster et al. [4] summarized that an ideal two-phase flow model and its numerical solution method should possess the following characteristics: (1) it should be robust in describing the morphological changes of complex phase interfaces; (2) it can calculate the surface tension accurately; (3) it can deal with the large density ratio and viscosity ratio problems; (4) it should be capable of handling multiphase flow at various spatial and temporal scales. However, at present, a single numerical simulation method for multiphase flow is unable to take all above requirements into account. Currently, the widely used numerical simulation approaches for multiphase flows can be classified into different types from different perspectives, as shown in

Table 1. Commonly used numerical simulation methods for multiphase flows.

No.	Classification	Representative Methods
1	interface modeling	interface-capturing method, interface-tracking method
2	fluid motion	Euler’s method, Lagrange method, hybrid method
.	…	…

. Among them, the interface modeling approach, including the interface capturing method and the interface tracking method, enjoys more popularity in engineering practice [5,6], which will be introduced in detail in this work.

According to the interface morphology, the interface-capturing method shown in Table 1 can be divided into sharp-interface method and diffuse-interface method. The main differences between them are qualitatively compared in Figure 1 [7]. Generally, the commonly used sharp-interface methods include volume of fluid (VOF) [8], level set (LS) [9,10], and others [11,12]. The fundamental theories and algorithms of VOF and LS methods are relatively mature, and they have been applied in many commercial software packages [13]. However, there exist several distinctions between VOF and LS methods regarding the calculation accuracy of phase parameters, such as the interface curvature and surface tension, the ability of smoothing when physical quantities at the interfacial region change violently, and the computational performance of mass conservation, etc. To remedy the shortcomings of these two methods, several hybrid approaches have been proposed in recent years, such as the coupled level-set and volume-of-fluid (CLSVOF) method [14], the hybrid particle level-set (PLS) method [15], and the volume-of-fluid and level-set (VOSET) method [16], etc. Despite the fact that these hybrid approaches can make up for the limitations of VOF and LS methods to a certain extent, such as ensuring the high accuracy of interface-curvature calculation and mass conservation characteristics, the programming process becomes more complex and cumbersome due to the increasing complexity. In addition, the program robustness is restricted, and the accurate calculation of surface tension is still not well solved.

It is worth noting that both VOF and LS methods belong to the sharp-interface method. In this type of method, the phase interface is assumed to have zero thickness and the physical properties of two phases change suddenly at the interface, which is obviously inconsistent with the real physical situation. In the 19th century, Rayleigh [17] and van der Waals [18] developed the diffuse-interface method by introducing the gradient theory of fluid interfaces based on the laws of thermodynamics. In the diffuse-interface method, the phase interface is considered to have a certain thickness and the physical properties of two bulk phases change continuously at the transition zone. On this basis, Korteweg [19] proposed the constitutive relation of interfacial stress tensor based on a free-energy functional of mass density or molar concentration. The phase-field model (PFM) is a kind of diffuse-interface method. Based on the diffuse-interface theory, PFM matches well with the real physical situation: the phase interface has a certain thickness and the physical properties of each bulk phase change continuously and rapidly in the transition region; the surface tension in traditional multiphase flow models equals to the body force distributed in the transition region. Compared with the sharp-interface method, PFM can offer many advantages from different points of view: (1) PFM is based on the energy variation and principle of least action, which has a complete physical background and allows the possibility to probe into the physical properties of the phase interface essentially. For example, PFM can deal with the problems of a moving contact line in multiphase flows, and the breakup and merging processes of a droplet with large deformations of interface morphology. (2) Although the derivation process of PFM is complex, the final expression of PFM is relatively simple and similar with that of the LS method. Thus, PFM is also called a physically motivated LS method or “conservative” LS method. It is easy to be extended from two-dimensional (2D) problems to three-dimensional (3D) problems. (3) The reconstruction of the phase interface by the volume fraction in the VOF method and the initialization of distance function in the LS method are not required in PFM. Therefore, PFM takes into account the merits of both VOF and LS methods. (4) PFM is capable of simulating the complex interface properties of fluids with different properties by defining the energy functional in the phase-field equation, such as polymer solutions, two-phase flows with soluble surfactant, etc. Therefore, it is of great significance to investigate the phase-field method.

In view of the advantages of the diffuse-interface method, PFM has emerged as a powerful and flexible tool in the numerical simulation of multiphase flows in recent years. A large number of researchers have devoted great effort to PFM from different aspects. In this paper, the main research progress and advances of PFM in the numerical simulation of multiphase flows in recent years are reviewed and analyzed, including the fundamental theory, mathematical models, numerical algorithms, and representative applications. The existing primary challenges and future research trends of PFM in the numerical simulation of multiphase flows are put forth and summarized.

2. Theory and Mathematical Models

In this section, the basic concepts, fundamental theory and mathematical models of PFM are introduced (see Figure 2).

2.1. Fundamental Theory

The fundamental theory of the PFM is based on the diffuse-interface theory of van der Waals [18]. It is considered that when two kinds of isothermal immiscible and incompressible fluids are in contact with each other, instead of being described by a geometric sharp interface, the two phases are separated by a thin layer of interface with a certain thickness. The physical properties of the two bulk phases are homogeneous, while the components of the two bulk phases are mixed with each other at the interface (transition zone) with a continuous and smooth change. The interfacial thickness is mainly dependent on the equilibrium between the molecule dynamics of two bulk phases, the value of which is approximately the same with molecules’ interaction range. The mixing effect between two phases and the diffuse interface are schematically shown in Figure 3 [20].

According to the diffuse-interface theory, a continuous scalar parameter $\varphi \left(x,t\right)$ is introduced in PFM to characterize the two phases and the interface, which is also called the order parameter,

$$\varphi \left(x,t\right)=\frac{{\rho }_1{C}_1-{\rho }_2{C}_2}{{\rho }_1{C}_1+{\rho }_2{C}_2}$$

(1)

where $\varphi \left(x,t\right)$ is the order parameter; ${\rho }_1$ and ${\rho }_2$ are the local density of bulk phases 1 and 2, while $C_1$ and $C_2$ are the corresponding local concentrations with $C_1+C_2=1$.

It can be seen from Equation (1) that the order parameter can characterize different fluid component regions. When the fluid is in the bulk phase region, the value of the order parameter is −1 or 1. When the fluid is in the interfacial region, the value of the order parameter is between −1 and 1. Thus, we have

$$\varphi \left(x,t\right)=\{\begin{array}{ll}1,& {\mathrm{C}}_2=0\begin{array}{c}(\mathrm{fluid}1)\end{array}\\ -1,& { \mathrm{C}}_1=0\begin{array}{c}(\mathrm{fluid}2)\end{array}\end{array}$$

(2)

According to the diffuse-interface theory, for a given phase field, the internal driving force of the phase distribution is the variation of the total free energy in the phase region. Cahn and Hilliard [21,22], respectively, derived the general equation describing the free energy of an intensity scalar field, i.e., the phase-field-governing equation, also known as the Cahn–Hilliard (CH) equation. Considering the non-uniformity of the phase field at the interface, Cahn modified the free energy density in the model to introduce the influence of spatial heterogeneity. The modified free energy density consists of two parts: (1) One is derived from the local material composition in the bulk phase. After expanding the free energy density formula by the Taylor series formula, only the first term, namely the main energy F(ϕ), is left in the free energy density formula when the local material components remain uniform. (2) The other is from the material heterogeneity $\nabla \varphi$ at the phase interface, which is also known as the ‘gradient energy’. As a result, the free energy density function of a phase field can be expressed as

$$w\propto \left(\varphi ,\nabla \varphi \right)$$

(3)

where $w$ is the free energy density of the phase field; ϕ is the order parameter; and $\nabla \varphi$ is the gradient energy.

The modified Helmholtz free energy function of the Ginzburg–Landau type can then be given by

$$\{\begin{array}{l}w\left(\varphi ,\nabla \varphi \right)=\lambda \left(\frac{1}{2}{\Vert \nabla \varphi \Vert }^2+F\left(\varphi \right)\right)\\ W\left(\varphi ,\nabla \varphi \right)={\displaystyle {\int }_{\Omega }\lambda \left(\frac{1}{2}{\Vert \nabla \varphi \Vert }^2+F\left(\varphi \right)\right)}\mathrm{d}x\end{array}$$

(4)

where w is the free energy density of the phase field; W is the total free energy of the phase field; $\lambda$ is the coefficient of mixed energy density; and F(ϕ) is the free energy of bulk phases. The first term on the right-hand side of Equation (4) is the free energy of the phase interface, which represents the nonlocal interaction between two phases, and its value denotes the phase affinity (mixing trend) changing with the concentration gradient of two phases. The second term on the right-hand side of Equation (4) is the free energy of the bulk phase, indicating the separation property of two phases, which acts as the main driving force of phase separation, and its value depends on the concentration of two phases. In the phase separation process, the total free energy of the phase field decreases with the development of fluid to a pure phase.

In the numerical simulation of multiphase flows, the double-well potential function shown in Figure 4 is usually selected as the free energy functional of bulk phases $F\left(\varphi \right)={\left({\varphi }^2-1\right)}^2/4{\epsilon }^2$ with ε being the capillary width of interface. As can be seen from Figure 4, F(ϕ) has a maximum value at ϕ = 0 and minimum value at ϕ = ±1, which is consistent with the energy distribution characteristics of bulk phases and phase interface.

Based on the diffuse-interface theory, the mixing between two-phase fluids is triggered by the gradient flow of chemical potential at the interface. That is, the chemical potential is the main driving force causing the phase separation and diffusion in the two-phase flow, and it plays a dominant role in the variation of free energy of the interface. According to the principle of variational calculus, the chemical potential of a heterogeneous multiphase flow system can be defined as the variation of free energy to the order parameter

$$\mu =\frac{\delta W}{\delta \varphi }=\frac{\partial w}{\partial \varphi }-\nabla \cdot \frac{\partial w}{\partial \nabla \varphi }$$

(5)

where μ is the chemical potential; W is the total free energy; w is the free energy density; and ϕ is the order parameter.

Then, the chemical potential can be further written as

$$\mu =\lambda \left\{\frac{1}{2}\frac{\partial \left({\Vert \nabla \varphi \Vert }^2\right)}{\partial \varphi }+\frac{\partial F\left(\varphi \right)}{\partial \varphi }-\nabla \cdot \left[\frac{1}{2}\frac{\partial \left({\Vert \nabla \varphi \Vert }^2\right)}{\partial \varphi }+\frac{\partial F\left(\varphi \right)}{\partial \nabla \varphi }\right]\right\}=\lambda \left(f\left(\varphi \right)-\Delta \varphi \right)$$

(6)

where $f\left(\varphi \right)=F^{\prime }\left(\varphi \right)=\left({\varphi }^3-\varphi \right)/{\epsilon }^2$.

When the interface is close to the limit of the sharp interface, based on f(ϕ), the relationship between the mixed energy density coefficient λ and the traditional surface tension coefficient σ can be determined, by which PFM can be related with the traditional sharp-interface method. This is also the mathematical expression of surface tension described by PFM as the volume stress distributed in the interfacial region with a certain thickness. Take the one-dimensional (1D) phase interface as an example, the relationship between λ and σ at the 1D phase interface can be derived as follows [23,24,25,26].

According to the phase field theory, the traditional surface energy equals the diffusive mixing energy in the region. Thus, the surface tension coefficient can be obtained by integrating the free energy in this region,

$$\sigma =\lambda {\displaystyle {\int }_{-\infty }^{+\infty }\left(\frac{1}{2}{\left(\frac{\mathrm{d}\varphi }{\mathrm{d}x}\right)}^2+F(\varphi )\right)}\mathrm{d}x$$

(7)

where $\sigma$ is the surface tension coefficient.

When the diffusive interface is at equilibrium, the chemical potential is zero.

$$\mu =\lambda \left(-\frac{{\mathrm{d}}^2\varphi }{\mathrm{d}{x}^2}+f(\varphi )\right)=0$$

(8)

It is easy to know that $\underset{x\to \pm \infty }{\lim }F(\varphi )=0$ and $\underset{x\to \pm \infty }{\lim }\frac{\mathrm{d}\varphi }{\mathrm{d}x}=0$. Multiplying $f(\varphi )=\frac{{\mathrm{d}}^2\varphi }{\mathrm{d}{x}^2}$ by $\frac{\mathrm{d}\varphi }{\mathrm{d}x}$ and integrating it over the interval (-∞, x), one can obtain

$${\displaystyle {\int }_{-\infty }^{x}f(\varphi )}\frac{\mathrm{d}\varphi }{\mathrm{d}x}dx={\displaystyle {\int }_{-\infty }^x\frac{{\mathrm{d}}^2\varphi }{\mathrm{d}{x}^2}}\frac{\mathrm{d}\varphi }{\mathrm{d}x}\mathrm{d}x$$

(9)

Then, we have

$$F(\varphi )={\displaystyle {\int }_{-1}^{\varphi }f(\varphi )}\mathrm{d}\varphi ={\displaystyle {\int }_{-\infty }^x\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{1}{2}{\left|\frac{\mathrm{d}\varphi }{\mathrm{d}x}\right|}^2\right)}\mathrm{d}x={\frac{1}{2}{\left(\frac{\mathrm{d}\varphi }{\mathrm{d}x}\right)}^2|}_{-\infty }^x=\frac{1}{2}{\left(\frac{\mathrm{d}\varphi }{\mathrm{d}x}\right)}^2$$

(10)

When the diffuse interface is at equilibrium, Equation (7) can be further derived based on Equation (10).

$$\begin{array}{ll}\sigma & =\lambda {\displaystyle {\int }_{-\infty }^{+\infty }\left(\frac{1}{2}{\left(\frac{\mathrm{d}\varphi }{\mathrm{d}x}\right)}^2+F(\varphi )\right)}\mathrm{d}x=\lambda {\displaystyle {\int }_{-\infty }^{+\infty }2\sqrt{\frac{1}{2}{\left(\frac{\mathrm{d}\varphi }{\mathrm{d}x}\right)}^{2}F(\varphi )}}\mathrm{d}x\\ & =\lambda {\displaystyle {\int }_{-1}^1\sqrt{2}}\sqrt{F(\varphi )}\mathrm{d}\varphi =\sqrt{2}\lambda {\displaystyle {\int }_{-1}^1\frac{1-{\varphi }^2}{2\epsilon }}\mathrm{d}\varphi \\ & =\frac{2\sqrt{2}}{3}\frac{\lambda }{\epsilon }\end{array}$$

(11)

where σ is the surface tension coefficient; λ is the mixed energy density coefficient; and ε is the capillary width of the interface.

It can be clearly seen from Equation (11) that when the capillary width of the interface ε approaches to 0, the mixed energy density coefficient λ must approach to 0 as well. The ratio of the two parameters determines the surface tension coefficient σ under the limiting condition close to the sharp interface. However, it is worth noting that Equation (11) is only valid at the equilibrium state.

2.2. Mathematical Models

In nature and engineering practice, multiphase flow often occurs under isothermal conditions; the influence of temperature change is generally ignored. Under isothermal conditions, the mathematical models of PFM include the phase-field equation, momentum equation and continuity equation, which are, respectively, introduced below.

2.2.1. Phase-Field Equation

Commonly seen forms of the phase-field equation include the conservative Cahn–Hilliard (CH) equation [21,22] and the non-conservative Allen–Cahn (AC) equation [27], both of which are derived by the principle of variational calculus. According to the second law of thermodynamics, the value of Helmholtz free energy is the smallest when the system is at an equilibrium state. Therefore, the variation of order parameter is determined by the minimization principle of Helmholtz free energy, and the phase-field equation can be expressed as

$${\varphi }_t+\left(u\cdot \nabla \right)\varphi =-M\frac{\delta W}{\delta \varphi }$$

(12)

where ϕ is the order parameter; u is the velocity; W is the total free energy; and M is the mobility parameter related to the relaxation time, which determines the intensity of the diffusion effect, that is, the mixing capacity of fluid molecules. In the two-phase flow, M determines the speed at which the phase interface reaches the equilibrium state.

The Cahn–Hilliard equation can be obtained by the variational operation to $\delta W/\delta \varphi$ in H⁻¹ space. It is a conservative phase-field equation that can preserve the conservation of order parameter in the whole computational domain for incompressible fluid. The expression of the Cahn–Hilliard equation is shown as below:

$${\varphi }_t+\left(u\cdot \nabla \right)\varphi =M\Delta \mu$$

(13)

If the variational operation is performed for $\delta W/\delta \varphi$ in L² space, the original formula of the Allen–Cahn equation can be obtained. It is a non-conservative phase-field equation given by

$${\varphi }_t+\left(u\cdot \nabla \right)\varphi =-M\mu$$

(14)

To keep the mass conservation in the numerical computation, a Lagrange operator is added to the Allen–Cahn equation. The modified Allen–Cahn equation can then be written as

$${\varphi }_t+\left(u\cdot \nabla \right)\varphi =-M\lambda \left(f\left(\varphi \right)-\Delta \varphi -\beta \left(\varphi \right)\right)\begin{array}{cc},& \beta \left(\varphi \right)=\frac{{\displaystyle {\int }_{\Omega }^{}f\left(\varphi \right)\mathrm{d}x}}{{\displaystyle {\int }_{\Omega }^{}\mathrm{d}x}}\end{array}$$

(15)

where the Lagrange operator $\beta \left(\varphi \right)$ can ensure $\frac{\mathrm{d}}{\mathrm{d}t}{\displaystyle {\int }_{\Omega }^{}\varphi \mathrm{d}x=0}$.

Except for above derivation, it is also possible to derive the PFM in a geometric manner starting from the LS method, in which the PFM is known as the so-called “conservative” LS method [28]. It should be noted that the above Cahn–Hilliard equation, Allen–Cahn equation and modified Allen–Cahn equation demonstrate different interfacial dynamics in the absence of external driving forces or advection [29,30]. Generally, the Allen–Cahn equation evolves by minimizing the mean curvature, the Cahn–Hilliard equation evolves through a mixture of bulk and surface diffusion depending on the mobility form, and the evolution of the modified globally conserved Allen–Cahn equation is somewhere between that of the Cahn–Hilliard and Allen–Cahn equations. Additionally, the choice of Cahn–Hilliard equation versus Allen–Cahn equation is dependent on a variety of factors such as the number of components in the fluid, whether there are additional driving forces between the phases beyond the interfacial energy, or whether there are reactions at the interfacial region, etc. Compared with the Allen–Cahn equation, the Cahn–Hilliard equation is more widely used in the numerical simulation of multiphase flows in wide applications. The primary challenges in solving the Cahn–Hilliard equation lie in the four-order term $\nabla$⁴$\varphi$ and the nonlinear term $f(\varphi )$.

2.2.2. Momentum Equation

The general form of momentum equation for fluid flow can be written as follows, whether it is a single-phase flow or a multiphase flow:

$$\rho \left(u_t+\left(u\cdot \nabla \right)u\right)=\nabla \cdot \left(\eta D\left(u\right)-pI+{\tau }_e\right)$$

(16)

where ρ is the fluid density; η is the dynamic viscosity of fluid; u is the velocity; $\eta D\left(u\right)$ stands for the viscous stress, $\eta D\left(u\right)=\eta \left(\nabla u+{\left(\nabla u\right)}^{\mathrm{T}}\right)$; $-pI$ denotes the pressure stress; and ${\tau }_e$ represents the external stress.

In the derivation of the phase-field equation, the phase interface is treated as a free condition. However, in the numerical calculation, it is essential to convert the surface tension into a body force for both the sharp-interface method and the diffuse-interface method. Therefore, a body-force term should be added to the momentum equation of PFM, the effect of which equals the surface tension. This treatment corresponds to the energy integral in the interfacial region. That is, PFM applies the interfacial free energy to describe the interfacial tension. Therefore, compared with traditional methods, PFM can better reflect the physical meaning of a problem.

(1)

Matched density;

When the density of two-phase fluids matches with each other, we have ${\tau }_e=-\lambda \left(\nabla \varphi \otimes \nabla \varphi \right)$. Substituting ${\tau }_e$ into Equation (16) yields the following momentum equation of PFM:

$$\rho \left(u_t+\left(u\cdot \nabla \right)u\right)=\nabla \cdot \left(\eta D\left(u\right)\right)-\nabla p-\varphi \nabla \mu$$

(17)

(2)

Non-matched density with a small density difference;

In the presence of non-match density, if the density difference between two phase fluids is small, the Boussinesq approximation can be used to treat the non-matched density as a uniform quantity. An additional gravity term is added in Equation (17) to introduce the effect of density difference. Then, the momentum equation of PFM can be expressed as

$${\rho }_0\left(u_t+\left(u\cdot \nabla \right)u\right)=\nabla \cdot \left(\eta D\left(u\right)\right)-\nabla p-\varphi \nabla \mu +g\left({\rho }_1,{\rho }_2\right)$$

(18)

where ${\rho }_1$ and ${\rho }_2$ are the local density of bulk phases 1 and 2; ${\rho }_0=\frac{{\rho }_1+{\rho }_2}{2}$; and $g\left({\rho }_1,{\rho }_2\right)$ is the additional gravity term introduced by the Boussinesq approximation, $g\left({\rho }_1,{\rho }_2\right)=\left(\frac{1-\varphi }{2}{\rho }_1+\frac{1+\varphi }{2}{\rho }_2\right)g$; g is the gravitational acceleration.

(3)

Non-matched density with a large density difference.

When the density difference of two phase fluids is large enough, the Boussinesq approximation is no longer applicable. Both the fluid density and viscosity are functions of the order parameter that need to be updated repeatedly in the calculation. Under this condition, the mass conservation cannot be guaranteed by the non-solenoidal velocity field ($\nabla \cdot u\ne 0$). To solve this problem, Guermond et al. [31] proposed that the two-phase fluid components can only be mixed by pure convection. The velocity at a certain point in the flow field is defined as the velocity of the phase component occupying this point. In this method, the momentum equation is modified physically to make the mixing of two phases driven only by convection, so as to guarantee the energy stability of the method. The momentum equation of PFM with a large density ratio can be written as follows based on the above pure convection mixing:

$${\sigma }_{\rho }{({\sigma }_{\rho }u)}_t+\frac{1}{2}\nabla \cdot (\rho u)u+\rho (u\cdot \nabla )u=\nabla \cdot \left(\eta D(u)\right)-\nabla p-\varphi \nabla \mu$$

(19)

where ${\sigma }_{\rho }=\sqrt{\rho }$; $\rho$ is the fluid density, $\rho =\frac{{\rho }_2-{\rho }_1}{2}\varphi +\frac{{\rho }_1+{\rho }_2}{2}$; and $\eta$ is the dynamic viscosity, $\eta =\frac{{\eta }_2-{\eta }_1}{2}\varphi +\frac{{\eta }_1+{\eta }_2}{2}$.

In addition, Abels et al. [32] developed a model that satisfies both the divergence of the velocity field and the energy law based on the idea of average volumetric velocity. In this case, a diffusion flux caused by the chemical potential and the density difference is introduced into the momentum equation. The modified momentum equation of PFM with a large density ratio is then given by

$$\rho \left(u_t+(u\cdot \nabla )u\right)+J_D\nabla u=\nabla \cdot \left(\eta D(u)\right)-\nabla p-\varphi \nabla \mu$$

(20)

where $J_D$ is the diffusion flux, $J_D=\frac{{\rho }_1-{\rho }_2}{2}M\nabla \mu$; $\rho$ is the fluid density, $\rho =\frac{{\rho }_2-{\rho }_1}{2}\varphi +\frac{{\rho }_1+{\rho }_2}{2}$; and $\eta$ is the dynamic viscosity, $\eta =\frac{{\eta }_2-{\eta }_1}{2}\varphi +\frac{{\eta }_1+{\eta }_2}{2}$.

2.2.3. Continuity Equation

Besides the phase-field equation and momentum equation, the multiphase flow should satisfy the continuity condition

$${\rho }_t+\nabla \cdot \left(\rho u\right)=0$$

(21)

It should be noticed that when the density of two phase fluids matches with each other or the density difference of two phase fluids is small, the fluid density can be approximately treated as a constant independent of time. However, when the density of two phases differs largely, the variation of density with time should be considered and the two phase fluids are treated as the compressible fluid.

2.2.4. Boundary Condition and Initial Condition

In addition to the phase-field equation, momentum equation and continuity equation, the corresponding boundary conditions and initial conditions should be given for a specific problem. Assuming that the computational domain is $\Omega$, the boundary conditions and initial conditions can be, respectively, presented by

$${u\cdot n|}_{\partial \Omega }=0,\hspace{1em}{\nabla \mu \cdot n|}_{\partial \Omega }=0,\hspace{1em}{\nabla \varphi \cdot n|}_{\partial \Omega }=0$$

(22)

$${u|}_{t=0}=u_0,\hspace{1em}{\varphi |}_{t=0}={\varphi }_0$$

(23)

where n denotes the unit normal vector of the boundary; $\partial \Omega$ is the boundary of the domain; and $u_0$ and ${\varphi }_0$ are the initial velocity and order parameter.

By coupling above phase-field equation, momentum equation and continuity equation, we can use the Cahn–Hilliard and Navier–Stokes (CH-NS) system or Allen–Cahn and Navier–Stokes (AC-NS) system to simulate multiphase flows. It should be noted that the coupled CH-NS is often referred as the Model H phase field in the broader applied physics communities. The classical Model H was initially developed by Hohenberg and Halperin [33] for simulating the binary incompressible two-phase flows with density-matched fluids, and was later generalized for simulating binary incompressible two-phase flows with variable density. Then, Lowengrub and Truskinovsky [34] extended Model H to be thermodynamically consistent for multiphase flows with large density ratios. And different modifications were also developed after Model H, such as the work of Guermond et al. [31], Abels et al. [32], and others [35,36].

3. Numerical Methods

In this section, the commonly used numerical schemes and algorithms used to solve the governing equations of PFM in the application of multiphase flows are presented (see Figure 5). Multiphase flow is a typical pressure–velocity coupling problem. Therefore, in the numerical simulation of multiphase flows by PFM, the coupling solution of the phase-field equation and momentum equation should be carried out to obtain the velocity, pressure, order parameter and surface tension, etc., with physical meanings. In the numerical calculation, the following should be noted: (1) Proper discrete schemes should be selected to discretize the coupled CH-NS system or AC-NS system. These discrete schemes should possess good numerical properties, such as the numerical stability, the conservation, etc. In addition, the effective treatment of the fourth-order term in the phase-field equation is one of the key issues in PFM. (2) An appropriate pressure–velocity coupling algorithm should be used to decouple the velocity and pressure in numerical simulations. The commonly used pressure–velocity coupling schemes are shown in

Table 2. Pressure–velocity coupling schemes.

Coupling Schemes	Representative Methods
All equations are solved simultaneously	Solving all variables simultaneously and globally;
	Solving part of variables simultaneously and globally;
	Solving all variables simultaneously and locally.
Equations are solved separately	Non-pressure-based approach;
	Vorticity–stream function method, vorticity–velocity method, etc.; Pressure-based approach: Pressure-correction method, projection method, fractional step method, artificial compression method, pressure Poisson equation method, etc.

[37,38]. Among them, the pressure-based approach is more widely used, such as the pressure-correction method, the projection method, and fractional step method, etc. In the application of these algorithms, the corresponding modifications are usually needed to obtain better numerical properties, such as good accuracy, robustness, mass conservation, and discrete energy decay, etc. One of the challenges of PFM is how to efficiently deal with the multiphase flow with large density and viscosity contrasts. (3) The selected algorithm should have a good applicability, which is capable of simulating the multiphase flows in both simple/regular and complex/irregular regions. (4) The numerical methods should be able to handle the multiphase flows at different spatial scales (such as the macroscale and microscale). With different points of concern, much research has been conducted on the numerical methods of PFM and some representative works are introduced and summarized below.

3.1. Discretization of Governing Equation

In the numerical solution of PFM, firstly, appropriate discrete schemes should be selected to discretize the phase-field equation and momentum equation. The discrete schemes should ensure the mass conservation, numerical stability and energy decaying of the discrete systems with time step as large as possible. In the field of computational fluid dynamics (CFD), a variety of discrete schemes for unsteady term, convection term, and diffusion term have been well developed. For example, the commonly used discrete schemes for unsteady term include the first-order forward difference scheme, second-order Adams–Bashforth scheme, second-order Crank–Nicolson scheme, and Runge–Kutta scheme, etc.; the widely used convective schemes contain the first-order upwind scheme, second-order upwind scheme, QUICK scheme, bounded scheme, and high-order compact scheme, etc.; the frequently used discrete schemes for the diffusion term include the second-order central difference scheme, the fourth-order central difference scheme, etc. [39,40].

It is worth noting that there is a fourth-order term of order parameter in the phase-field equation. Generally speaking, most of convection–diffusion equations are second-order partial differential equations (PDEs); thus, the treatment of this fourth-order term is the focus in the discretization of phase-field equation. If the fourth-order term is discretized directly, the order parameter at one grid point to be solved will be related to that at multiple surrounding grid points, which not only seriously affect the calculation efficiency, but also deteriorate the numerical stability. The time step will also be restricted to a small value. Therefore, it is necessary to apply some special treatments to deal with the fourth-order term carefully, among which the convex-splitting technique and introducing a stabilizing term are commonly used [39,40]. Especially, the convex-splitting technique is popular to reformulate it as a pair of coupled second-order equations: one for the chemical potential, and the other for the evolution of order parameter. For example, suppose the free energy function of bulk phase is a double-well potential function, i.e., $F\left(\varphi \right)={\left({\varphi }^2-1\right)}^2/4{\epsilon }^2$, then $F\left(\varphi \right)$ can be split into two different convex functions $F_c\left(\varphi \right)$ and $F_e\left(\varphi \right)$ by the convex-splitting method,

$$F_{\mathrm{c}}\left(\varphi \right)=\frac{4{\varphi }^2+1}{4{\epsilon }^2},\hspace{1em}{F}_{\mathrm{e}}\left(\varphi \right)=\frac{6{\varphi }^2-{\varphi }^4}{4{\epsilon }^2}$$

(24)

where $F\left(\varphi \right)=F_{\mathrm{c}}\left(\varphi \right)-F_e\left(\varphi \right)$, $F_{\mathrm{c}}^{″}\left(\varphi \right)\ge 0$ and $F_{\mathrm{e}}^{″}\left(\varphi \right)\ge 0$. By using the convex-splitting scheme, the resulting discrete system of Cahn–Hilliard and Allen–Cahn equations is still unconditionally stable and satisfies the discrete energy law. Readers of interest can refer to [7,41,42] for the detailed proof.

Based on the convex-splitting method, many discrete schemes with good numerical properties have been developed for the discretization of PFM. For example, Shin et al. [40] combined the convex-splitting method with a multi-stage implicit–explicit Runge–Kutta method, and developed a convex-splitting Runge–Kutta (CSRK) scheme with high-order accuracy in time. The CSRK scheme can solve the Allen–Cahn, Cahn–Hilliard and phase-field crystal equations well with the calculation results possessing three-order accuracy. By using the adaptive mesh refinement and a nonlinear multigrid finite-difference method, Chen et al. [43] developed efficient energy stable numerical methods for solving the Cahn–Hilliard system with isotropy and strong anisotropy based on the convex-splitting approach. Liu et al. [44] proposed a fully discrete virtual element scheme coupling the virtual element discrete nonlinear term with the convex-splitting method for the mixed form of Cahn–Hilliard equation in arbitrary polygons. The proposed scheme was proved to possess properties such as the unconditional unique solvability, mass conservation, discrete energy decaying and large time step. Bu et al. [45] presented a stable second-order numerical scheme for space-fractional Cahn–Hilliard and Allen–Cahn equations, which performed the spatial discretization through the Fourier spectral method and the temporal discretization by the convex-splitting method, respectively. Numerical experiments showed that the scheme can preserve mass and is unconditionally energy stable.

In addition to the convex-splitting method, some other accurate and efficient numerical schemes have been put forward for the phase-field equation. For instance, Aboelenen et al. [46] developed a high-order nodal discontinuous Galerkin method for the linearized fractional Cahn–Hilliard equation. Through this method, the derivatives of fractional order in space can be expressed as a composite of first-order derivatives and integrals, and the linearized fractional Cahn–Hilliard system can be transformed as a low-order system of differential or integral equations. Zhao et al. [47] proposed a second-order fully discrete linear scheme for a thermodynamically consistent hydrodynamic PFM of binary compressible viscous fluids, in which the Crank–Nicolson scheme was firstly applied to discretize the model in time and a semi-discrete partial differential equation (PDE) was output, and then the second-order finite-difference scheme was used to discretize the semi-discrete PDE on staggered grids to obtain a fully discrete scheme in space. Yang et al. [48] put forward a second-order time-marching scheme based on an invariant energy quadratization method for the viscous Cahn–Hilliard equation with hyperbolic relaxation. In this method, all nonlinear terms were discretized by the semi-explicit scheme to arrive at a linear system of a symmetric positive definite equation. Numerical cases illustrate the proposed scheme is unconditionally energy stable and can converge in time with the second-order accuracy. Ahammad et al. [49] and Alam [50] presented a wavelet-based approach for the phase-field modeling of two-phase flows, in which the interfacial dynamics are described by the AC-NS equation and are solved by a weighted residual collocation method based on Deslauriers–Dubuc interpolating wavelets.

3.2. Multiphase Flows with Large Density and Viscosity Contrasts

In actual multiphase flows, it is common to see that the physical properties of two phase fluids differ greatly. For example, the densities of the oil phase and gas phase can up to thousands of orders different, which poses a great challenge to the stability of numerical algorithms. Therefore, another difficulty in using PFM to simulate multiphase flows is to deal with the fluid phase with large density and viscosity contrasts. The difference exists in numerical algorithms for handling multiphase flows with matched density and variable density mainly lies in the momentum equation. Due to the influence of variable density and viscosity, the numerical algorithm for non-matched-density multiphase flows is usually difficult to preserve the properties of divergence-free velocity field, mass conservation and energy stability, etc.

At present, there are two commonly used numerical methods for non-matched-density flows with large density and viscosity ratios. One is to modify the momentum equation based on a pure convection mixing idea, in which the fluid mixing is only driven by the convection and then the energy stability of numerical algorithms can be ensured. For example, Ding et al. [51] applied the above algorithm to study the applicability of PFM for incompressible two-phase flows with large density and viscosity differences; the convergence properties and numerical accuracy of the algorithm for binary mixtures with large density ratio was investigated. Kay et al. [52] proposed an efficient multigrid finite-element solver for the variable density CH-NS system based on this algorithm, in which the nonlinear multigrid method was used to solve the Cahn–Hilliard equation, and the multigrid preconditioned GMRES approach was used to solve the Navier–Stokes equation. Other cases can also be found in the literature [53,54,55], etc. The other way of dealing with multiphase flows with large density and viscosity contrasts is to treat the mixed fluids as a compressible fluid, and the interfacial velocity is calculated by the volume average or the mass average. The commonly used method is based on an additional flux caused by the chemical potential and density difference. For instance, Dong [56] developed an efficient numerical method in the CH-NS system for wall-bounded flows of immiscible incompressible two-phase fluids with a large density ratio, in which the dynamic contact angle boundary conditions were applied. Gao et al. [57] and Yu et al. [58] presented energy-stable numerical schemes for a moving contact line problem of two-phase flow with variable densities and viscosities; the coupled system of Cahn–Hilliard equation with dynamic contact line conditions and Navier–Stokes equation with generalized Navier boundary conditions was solved efficiently. Other representative works can be found in the literature [36,59,60,61].

In addition to the two widely used methods introduced above, several other numerical approaches have been developed to overcome the challenges in solving multiphase flows with large density and viscosity ratios. For example, Joshi and Jaiman [61] developed a positivity-preserving and conservative variational scheme for the phase-field modeling of two-phase flows with large density differences. The mass conservation was achieved by a Lagrange multiplier and the energy stability was guaranteed by the mid-point approximation method. Wang et al. [62] proposed an entropy–viscosity method (EVM) for the Cahn–Hilliard system to simulate two-phase flows with large density and viscosity ratios at a high Reynolds number. Through the artificial interface-compression method (AICM), which is based on the EVM, the stability of the simulation process was realized and the sharpness of the phase interface was well maintained. The readers can be referred to [63,64,65] for other representative works.

3.3. Multiphase Flows in Complex/Irregular Domains

Researches on PFM in the early stage mainly focused on the basic theory and numerical algorithms, most of which were aimed at the problems in simple or regular geometries. Few studies were conducted for the problems in complex or irregular domains. In recent years, with the development of phase-field theory and the advances in numerical algorithms, PFM has been frequently applied to complex domains with irregular boundary shapes [66,67,68,69,70,71], in which the description of domain geometry, mesh generation and algorithm stability are the primary challenges that attract much attention.

For example, Li et al. [66] developed a PFM applicable to the multicomponent Cahn–Hilliard system in complex domains. With the aid of an adaptive mesh refinement technique and an unconditionally gradient-stable nonlinear splitting scheme, the robustness of the proposed model was guaranteed. Luo et al. [67] solved the CH-NS model with generalized Navier boundary conditions by using the finite-element method (FEM) based on unstructured mesh, and 3D droplet spreading on a rough solid surface was investigated. In this method, the mass compensation algorithm was applied to guarantee the mass conservation of droplets and the domain decomposition approach was used to improve the computational efficiency. Li et al. [68] proposed a simple direct discretization method for the Cahn–Hilliard equation on an evolving surface, in which the unstructured triangular mesh was used to discretize the curved surface and the surface evolution was realized by moving the mesh nodes according to the given velocity field. Numerical experiments illustrated that the proposed method has second-order accuracy and is easy to implement. This method was also applied to solve the phase-field crystal equation on surfaces [69]. Zimmermann et al. [70] put forward an iso-geometric FEM for mass-conserving phase transitions on the deforming surface. The Cahn–Hilliard equation was used to describe the phase transition on the curved surface, and the Kirchhoff–Love thin-shell equation was applied to capture the surface deformation. In addition, the surface geometry was described by the structured NURBS meshes and unstructured quadrilateral meshes. Guo et al. [71] coupled a diffuse domain (DD) method with a quasi-incompressible CH-NS model to simulate two-phase flows with moving contact lines in complex geometries. In the DD method, the original complex domain is extended to a larger regular domain, and the complex domain boundary is replaced by an interfacial region with finite thickness. By using the DD method, no modifications are required in standard commercial software packages and the stationary/moving complex boundary can be imposed easily.

3.4. Multiphase Flow at Microscale

Besides the application in simulating multiphase flows at macroscale, PFM is also a powerful tool for multiphase flow simulations at mesoscale when combined with the lattice Boltzmann method (LBM). The primary advantages of PFM within the framework of LBM can be summarized as follows: (1) The mass and heat transfer of multiphase flows can be described by LBM well without the restriction of a continuous hypothesis. (2) The evolution object of LBM is the distribution function of discrete particles. Thus, LBM can not only easily describe the interactions between different fluids or between the fluid and solid, but also can conveniently deal with various complex boundaries. (3) The collision process of discrete particles in LBM occurs partially in the domain, which is in general more efficient and easier to be parallelized.

In recent years, many scholars have focused on the phase-field–lattice Boltzmann method (PF-LBM), in which different distribution functions were used to describe the phase field, velocity and pressure [72,73,74,75,76,77,78,79,80,81,82,83,84,85,86]. By coupling LBM with PFM, it can be applied to solve multiphase flows with large differences in physical properties, and various measures are taken to preserve the mass conservation of different phases. For example, Fakhari et al. [72] presented a multiple-relaxation-time LBM (MRT-LBM) for low-viscosity binary fluids, which was then coupled with PFM for multiphase flows with moderate density ratios. The case study showed that PF-LBM can easily handle the phase interface with a large deformation rate. Banari et al. [73] developed an efficient PF-LBM for 3D multiphase flows based on the Cahn–Hilliard system. Both the phase-field equation and the momentum equation were solved by LBM in the framework of a parallel GPGPU co-processor for acceleration. This method can solve two-component multiphase flows with large density ratios at a high Reynolds number efficiently and accurately. Liang [74] established a PF-LBM for 2D and 3D multiphase flows in complex microchannels. The presented model can achieve an overall improvement in the accuracy and stability of the capturing interface, and the macroscopic pressure and velocity can be calculated much easier. Geier et al. [75] proposed a single-relaxation-time LBM (SRT-LBM) based on an improved Allen–Cahn model, which can recover the conservative phase-field and guarantee local and global mass conservation. Wang et al. [76] developed an MRT-LBM for a 3D multiphase flow with moderate density ratios, and the method adopted two distribution functions: one is the order parameter distribution function for solving the Cahn–Hilliard equation, and the other is the pressure distribution function for solving the momentum equation. Fakhari et al. [77] proposed a mass-conserving MRT-LBM for immiscible two-phase flows, in which the interface tracking and pressure evolution were described by the PF-LBE. This method has also been used to model 3D contact line dynamics on curved boundaries [78], and is capable of simulating multiphase flows at high density ratios and high Reynolds numbers in the Cartesian coordinate. Su et al. [79] applied an improved Allen–Cahn PF-LBM to numerically study the rising dynamics of a single bubble with large density ratios in quiescent viscous fluid, and the influences of Eotvos/Reynolds number, density ratio, viscosity ratio, bubble size and bubble shape were systematically investigated. Zhang et al. [80] conducted a comparative study on the capability of PF-LBM and the Lagrange–Euler finite-element method (ALE-FEM) in solving multiphase flows. Numerical experiments demonstrated that the PF-LBM is more dissipative than the ALE-FEM due to the nature of PFM. Zhang et al. [81] developed an LBM based on the PFM for simulating large-density-ratio two-phase flows. An improved MRT-LBM was developed to solve the conserved AC equation. By utilizing a nondiagonal relaxation matrix and modifying the equilibrium distribution function and discrete source term, the conserved AC equation can be correctly recovered by the proposed MRT-LBM with no deviation term. Other works can also be referred to [82,83,84,85].

It is worth pointing out that the influence of temperature change is usually neglected in most of the literature when the PFM is used to study multiphase flows. If the effect of temperature variation is taken into account, it is necessary to consider the influence of temperature on the physical properties of phase fluids, surface tension, etc. However, the relevant research is still lacking at present. Hu et al. [86] proposed a PF-LBM based on the conservative Allen–Cahn equation for thermocapillary flows. Three different distribution functions were applied to describe the evolution of order parameter, velocity and temperature, respectively. The proposed PF-LBM possesses high numerical stability and is capable of simulating thermocapillary flows with large density ratios and thermophysical properties contrasts such as specific heat capacity, thermal conductivity, etc. As the standard LBM is only applicable to Cartesian grids, the vast majority of research on PF-LBM focuses on multiphase flows in Cartesian coordinates, and that for complex geometries are rarely studied. Ambrus et al. [87] put forward a finite-difference PF-LBM based on the Cahn–Hilliard system, which was able to solve multicomponent multiphase flows on curved surfaces by rewriting the Boltzmann equation through vielbein formalism on arbitrary geometries. This method can be easily extended to other curved surfaces and coupled to other dynamic phenomena. In addition to the traditional LBM, some scholars have proposed energy-based MRT-LBM. For example, Alpak et al. [88] presented a PF-LBM based on the Helmholtz free energy minimization that couples the Cahn–Hillard equation and the energy-based LB equation. By this model, the DNS of two-phase flows at a pore-scale is achievable on 3D digital images of real rocks obtained by micro-CT scanning.

3.5. Multi-Component Multiphase Flows

The PFMs are also capable of solving the multi-component multiphase flow problems. The solution of this type of problem is more complex than multiphase flows because of the following reasons: (a) the relations between phases and components should be considered; (b) the component equations should be modeled; and (3) the efficient numerical schemes and algorithms are more complex and challenging to develop. Some scholars are devoted to this topic. For the state of the art of the of PFMs for multi-component systems, Nestler et al. [89] reviewed the application of PFM for studying phase transitions and microstructure evolution in multi-component systems. Kim [90] presented a review on PFMs and their numerical methods for multi-component fluid flows with interfacial phenomena. Santra et al. [91] comprehensively analyzed the development of PFM for multicomponent and multiphase flows in microfluidic systems. For the research work, Shi and Wang [92] put forth a PFM consisting of CH-NS equations for the simulation of three-component immiscible flows on a solid surface. The core idea of the proposed approach was to develop a particular formulation of boundary conditions to satisfy the extra consistency conditions for the system. Huang and Lin [93] developed a consistent and conservative model for multiphase and multicomponent (namely N-phase-M-component) incompressible flows based on the PFM. In this model, components were dissolvable in some specific phases and a dissolvability matrix was defined to indicate relations between phases and components. The component equation was finalized by incorporating the thermodynamical effect from the PFM so that the volume fraction equation derived from the component equation was consistent with the one from the PFM. Far et al. [94] proposed a PFM to simulate the soluble surfactant in a multicomponent multiphase system. In their model, the LBM was used to solve the Cahn–Hilliard equation for the order parameter and the surfactant concentration, and the interfacial tension of the system in the presence of the surfactant can be calculated analytically. Hester et al. [95] developed a second-order PFM to combine melting and dissolution in multi-component flows. In this method, a general framework for the asymptotic analysis of diffuse-interface methods in arbitrary geometries was put forth to eliminate the error proportional to the thickness of the phase interface. Yang [96] proposed PFMs for the evaporation of a single droplet and multiple droplets on a substrate. In the proposed model, to describe the evaporation of multiple droplets under different contact angles or pining boundary conditions, the multi-component Allen–Cahn model was adopted in the computation.

3.6. Others

The discretization and solution of PFM can be realized by traditional numerical methods such as the finite-difference method, finite-element method, and finite-volume method, etc. Other effective numerical methods have also been developed by considering the specific characteristics of the investigated problems. For example, Ahammad et al. [49] and Alam [50] solved two-phase flows described by the AC-NS equation through the wavelet transform method.

In PFMs, the interface width and the spatial/temporal steps should be chosen carefully and their values are always small enough in order to obtain accurate simulation results. Therefore, limited by computing capabilities and the problem scale, at present, the applications of PFM are more common in mechanism studies than in real engineering practice. Barrett et al. [97] compared the practicality and accuracy of PFMs and sharp-interface models through Stefan problems with applications to crystal growth. It was found that sharp-interface models were computationally cheaper than PFMs for simulating interface evolutions in materials science. Egger et al. [98] and Cervera et al. [99], respectively, compared the PFM with the family of FEM methods (XFEM, GFEM, SBFEM, MFEM) in the application of linear elastic fracture mechanics and quasi-brittle cracking. Bazant et al. [100] and Diehl et al. [101] compared the PFMs with the peridynamic model and crack band model in the study of fracture tests and engineering fracture mechanics. Jain [102] proposed an accurate conservative PFM for the simulation of two-phase flows; it was showed that the proposed model imposes a lesser Courant–Friedrich–Lewy restriction and is thus less computationally expensive compared to previous conservative PFMs. Overall, the cost-efficiency issue of the PFM is challenging.

On the premise of ensuring the numerical accuracy and energy stability of the algorithm, many scholars have focused on highly efficient solvers or acceleration techniques for solving the discrete PFM system, such as the adaptive mesh and time step, multigrid method, domain-decomposition technique, and CPU/GPU parallel computing, etc. Especially, with the improvement of GPU computing power and the GPU parallel programming framework, the multi-GPU parallel computing on high-performance clusters has become a significant means for the highly efficient simulation of PFM. For example, Zhao et al. [103] developed a first-order energy-stable efficient scheme for simulating binary mixtures of viscous fluids and nematic liquid crystals in the framework of PFM. The 3D efficient simulation was achieved by the CUDA on GPUs. Liu [104] carried out the CPU+GPU heterogeneous parallel computing of the phase field, flow field and other external physical fields when solving the multi-field coupled-phase PF-LBM, as shown in Figure 6. Xu [105] explored the parallel computing of PFM based on an adaptive FEM, and realized the phase-field simulation of Al-Cu binary alloys by using the parallel programming of MPI and MPI+OpenMP on a computer cluster. Alpak et al. [88] adopted the GPGPU parallel computing to efficiently simulate two-phase pore-scale flows described by the Helmholtz free energy—minimizing phase-field model within the framework of LBM. In general, the CPU/GPU parallel computing is an effective and universal acceleration method for large-scale simulation of multiphase flows in the framework of PFM.

It should be noted that, actually, the phenomena and applications of multiphase flow systems are complex. There are also other simulation methods that are widely used, such as the family of discrete particle methods [106,107,108] for dense particle systems, the methods with point particle approximation [109,110] for dispersed multiphase flows, etc. In this paper, we only mainly focused on the PFM. Through the above six subsections, the numerical methods for PFM were reviewed and analyzed from different perspectives. Although PFM has some obvious advantages, it should also be noted that there are some well-known pitfalls of the PFM in the numerical methods, such as the mesh anisotropy, artificial numerical trapping and pinning of the interface when it is insufficiently resolved or when there are large discontinuities in the material properties.

4. Applications

In recent years, PFM has been widely applied in a wide range of fields, such as hydrodynamics, petroleum science, materials science, and biomedicine, etc. Specifically, PFM is a valuable tool for the interface tracking of multiphase flows with large density or/and viscosity contrasts, and for the description of microstructure evolution of metals or alloys. Several typical applications of PFM in fluid mechanics and materials science are briefly reviewed here, especially in multiphase flows.

4.1. Application of PFM in Fluid Mechanics

In fluid mechanics, the primary difficulty of multiphase flow simulations is to accurately capture the evolution of the phase interface and calculate the surface tension, while the PFM possesses advantages in dealing with the above key issues. Nochetto et al. [111] established a mathematical model to describe flow behaviors of two-phase ferrofluid flows based on PFM, and presented an energy-stable scheme for the efficient solution of this model. Li et al. [112] developed a two-phase flow model considering the dynamic wetting effect based on PFM and a modified dynamic contact-angle model. The dynamic process of droplets impacting on a wall was simulated using this model on the OpenFOAM platform. Feng et al. [113] simulated oil–water two-phase flows in pore channels of porous media using the PFM. The effects of displacement velocity, fluid properties, and wettability on the oil recovery were analyzed in detail. In this study, the water flooding process of water–wet-sandstone porous media obtained from the proposed model was presented.

Balashov et al. [114] investigated the compressible isothermal two-phase flow of binary mixtures with surface effects based on a 1D CH-NS regularized system, and a semi-discrete scheme satisfying the thermodynamical consistency was proposed. Taking into account the double-sided model with equal diffusivity in liquid and solid phases, Subhedar et al. [115] extended the thin-interface limit of PFM to transport through melt convection. Conti et al. [116] discussed the well-posedness of PFM for the motion of binary fluids with large viscosity contrasts. In a 3D-bounded domain, the existence and uniqueness of the global weak solution and strong solution were proved for the Brinkman–Cahn–Hilliard system with logarithmic free energy density. Soligo et al. [117] combined the direct numerical simulation (DNS) of turbulent flow and the second-order-parameter PFM to depict the dynamic behaviors of droplets and surfactants. The breakage and coalescence events of surfactant-laden droplets in turbulent flows were simulated, and the influences of surface tension and surfactant types on these phenomena were analyzed. Zhang et al. [118] established a gas–liquid two-phase model by PFM to study the coalescence and separation of bubbles in an electrochemical system, and the coalescence of complex microbubbles on and departing the electrode and the separation of bubbles from the electrode surface were predicted. Based on a thermodynamically consistent CH-NS model, Bonart et al. [119] employed the PFM to simulate the commonly seen sliding phenomenon of liquid droplets along the solid surface, such as raindrops on leaves. Zhang et al. [120] studied the spontaneous behavior of droplets on partially wetting surfaces, and the critical radius below which the droplet will eventually disappear was theoretically predicted and numerically verified with different contact angles, domain sizes, and interface widths. Yang et al. [121] developed a novel PFM to describe compound droplets in contact with a solid, in which the solid was represented by the fixed phase of a four-component CH system. Comparison validated that the numerical and analytical results for the compound droplets on a flat solid are in good agreement with each other.

4.2. Application of PFM in Material Science

In material science, PFM can be conveniently applied to describe the crack propagation of materials [122,123,124,125,126,127,128] and the evolution of metal/alloy microstructures [129,130,131,132,133,134]. For example, Zhou et al. [122] numerically studied the fluid-driven dynamic crack propagation in poroelastic porous media by PFM, in which the arbitrary crack was controlled by the evolution of the phase field. The 2D and 3D examples of dynamic crack branching and their interaction with pre-existing natural fractures were carried out by the COMSOL Multi-physics software. Xia et al. [123] employed an extended PFM to simulate the hydraulic fracturing and cracking in heterogeneous fluid-saturated porous media. The evolution of crack phase fields during the hydraulic fracturing of realistic heterogeneous medium was presented. In Li et al. [124], PFM was used to simulate the crack nucleation and propagation in rock-like materials with prefabricated cracks. The comparison between simulation results and experimental data illustrated that PFM is capable of predicting the crack nucleation and propagation with different rock bridge angles under uniaxial compression. Ambati et al. [125] developed a novel PFM for modeling the ductile fracture of elasto-plastic material in a quasi-static linear regime. The merit of this model is that it couples the evolution of crack phase field with the accumulation of plastic strains in a thermodynamically consistent way. Makvandi et al. [126] proposed two different phase-field fracture models based on the strain gradient elasticity theory. A second-order and a fourth-order phase field formula were used, respectively, which eliminated the singularities at the crack tips and improved the performance of conventional models. Goswami et al. [127] presented a new adaptive PFM based on the dual mesh to deal with the brittle fracture, in which the elastic field adopts the coarser mesh and the phase field uses the finer mesh, respectively. To facilitate the information exchange between the meshes, an efficient data transmission algorithm was developed. Zhu et al. [128] put forward a PFM for the propagation of electrical tree in nanocomposite materials based on an improved Allen–Cahn equation. Simulation results demonstrated that the PFM can vividly reproduce the dendritic shape of electrical tree.

Wang and Cai [129] adopted the PFM to investigate rising behaviors of a single Taylor bubble in liquid lead–bismuth eutectic (LBE) along a vertical tube. This study proved that the PFM provides a useful alternative for simulating two-phase flows of Taylor bubble rising in the liquid metal. Using the PFM of a Karma thin interface with solid diffusivity, Zhang [130] employed a 3D adaptive mesh method to explore the influences of interface width and thermal noise on the dendrite growth in undercooled melt of Al-Si alloy. Pinomaa et al. [131] proposed a numerical approach that couples the CH-NS system with the heat transfer equation for modeling the powder-bed scale melt pool, in which the Cahn–Hilliard model was applied to predict the formation of rapid solidification microstructures for Ti-6-4 alloy in epitaxial and equiaxed growth regimes. The phase separation in liquid crystal was simulated in the work of Riahinasab et al. [132] using the Cahn–Hilliard equation, and the effect of cooling rate on the partitioning of nanoparticles into isolated domains was explored qualitatively. Zhang et al. [133] developed a novel phase-field crystal model based on the L²-gradient flow and put forward a second-order unconditionally energy-stable scheme with the mass conservation for solving this model. By comparing with the phase field of classical H⁻¹ gradient flows, the robustness of the new model was proved and the stability of the newly developed scheme was verified. Restricted to the article length, readers of interest can also refer to reference [134] for other representative works.

4.3. Other Interesting Applications

The PFM also enjoys popularity in many other scientific fields and engineering applications. For example, PFM can be used for binary image inpainting based on the modified Cahn–Hilliard model [135], and to predict the hydrate growth [136] and cell growth [137,138,139], and to model the formation of biofilm structure [140], etc. It should be mentioned that the application of PFM is not only restricted to small-scale structures; instead, it can also model the formation and evolution of large-scale structures and objects. For instance, the irregular structure of Saturn’s B ring was modeled by PFM based on the Cahn–Hilliard equation in the literature [141].

5. Conclusions and Future Work

Compared with conventional sharp-interface methods, the PFM possesses obvious advantages in the numerical modeling and simulation of multiphase flows in a wide range of applications. At present, the fundamental theory, mathematical model and numerical methods of PFM have been well studied by scholars at home and abroad. The PFM is more and more widely applied in scientific and engineering fields involving fluid mechanics, material sciences, etc. In this paper, the state-of-the-art of PFM in the numerical simulation of multiphase flows is briefly reviewed.

(1)

PFM is relatively mature in either fundamental theories or mathematical models. By coupling the commonly used phase-field models (Cahn–Hilliard or Allen–Cahn equations) or its modified versions with the momentum equation, the multiphase system of fluid/solid mixtures with large density and viscosity contrasts can be well modeled. However, it is also worth noting that classical PFM is always restricted to isothermal conditions without considering the influence of temperature variation, which needs more to pay attention to more effects in future study.

(2)

In terms of numerical methods, the discretization of the fourth-order term of order parameter and the numerical algorithm of multiphase flows with large density and viscosity ratios have been well developed. At present, the PFM is capable of dealing with multiphase flow problems with large differences in physical properties. However, the primary challenge of PFM in the numerical simulation is still to develop the robust, high-accurate and energy-stable numerical algorithms or schemes with mass conservation and thermodynamic consistency, especially to be applicable to the conditions of large density and viscosity ratios at a high Reynolds number. In addition, due to the limitation of computing capability, PFM now is most widely used in mechanism investigations, and its application in engineering practice is still at an initial stage. Thus, a stable and efficient numerical method for PFM is still a research hotspot for future work.

(3)

With respect to the application, PFM has been widely applied in various scientific and engineering fields such as fluid mechanics, material science, computer science, petroleum engineering, chemical engineering, biomedicine and astrophysics, etc. However, the application of PFM primarily remains on the mechanism study on simple or/and regular geometries; applications in practical engineering problems in complex/irregular domains are in initial stage and still need more effort. Thus, extending PFM to the real complex or large-scale engineering problems is still a research trend in the near future.