Unlike the phase-field model of the two-phase flow system that requires only one phase-field variable to represent the volume fraction of two phases, the Cahn-Hilliard phase-field model of three-phase incompressible flows usually requires three phase-field variables to formally represent the volume (or mass) of each component [1][2][3][4]. In the total free energy, the hydrophilic-hydrophobic tendency of each phase-field variable is independent, but to ensure the socalled free-leakage condition of the system, a Lagrange multiplier needs to be added in each Cahn-Hilliard equation, thereby establishing a highly nonlinear and coupled system. Not only that, for some specific physical phenomena, such as the so-called "total spreading" situation, but it is also necessary to add some coupled higher-order terms to the total free energy to ensure the well-posedness of the entire system. Finally, after combining the Navier-Stokes equation describing the characteristics of the incompressible flow field, a highly coupled and nonlinear complex dynamical system is obtained. Among them, the coupling terms can be divided into two categories, one is formulated introduce the numerical scheme, explain its implementations in detail, and prove its solvability and discrete energy dissipation law rigorously. In Section 4, many numerical examples are given to illustrate the accuracy and efficiency of the scheme. We finally give some concluding remarks in Section 5.
We briefly outline the Cahn-Hilliard phase-field system, which simulates the three immiscible fluid components proposed in [2][3][4]. A domain Ω is an open bounded, connected, subset in R d of d = 2, 3, with a sufficiently smooth boundary. We assume φ i (i = 1, 2, 3) to be the ith phase-field variable which represents the volume of the ith component in the fluid mixture, such that
1 inside the ith component, 0 outside the ith component, (
where x ∈ Ω , t is time in [0, T ]. A smooth layer with the thickness ϵ is used to connect the interface between 0 and 1. Assuming the mixture being perfect (called as free-leakage condition), the three unknowns φ 1 , φ 2 , φ 3 are linked by the following relationship:
There are several ways to promote the two-phase model to the three-phase case, cf. [1][2][3][4]. In this article, we use the model developed in [2][3][4] and assume the total free energy as
where ϵ is the order parameter to characterize the interfacial width, L(φ 1 , φ 2 , φ 3 ) is the linear part, and F(φ 1 , φ 2 , φ 3 ) is the nonlinear part. The linear part is given as
where the coefficient Σ i represents the "spreading" coefficient of the ith at the interface between jth phase and kth phase. One can deduce that the following conditions hold between the three surface tension parameters σ i j (σ 12 , σ 13 , σ 23 ) and three spreading coefficients (Σ 1 , Σ 2 , Σ 3 ) (see [2][3][4]):
Σ i = σ i j + σ ik -σ jk , i = 1, 2, 3.
(2.5)
Note Σ i might not be always positive. If Σ i > 0, it means that the spreading is "partial", if Σ i < 0, it is called "total". The nonlinear potential F(φ 1 , φ 2 , φ 3 ) is given as:
where Λ is a non-negative constant. Since φ 1 , φ 2 , φ 3 satisfy the free-leakage condition (2.2), F(φ 1 , φ 2 , φ 3 ) can be rewritten as
where
(2.8)
Regarding the total free energy and surface tension coefficients, some remarks are given as follows (see the details in [4]):
Remark 2.1. Let σ 12 , σ 13 and σ 23 be three positive numbers and Σ 1 , Σ 2 and Σ 3 defined by (2.5). For any Λ > 0, the bulk free energy F(φ 1 , φ 2 , φ 3 ) defined in (2.7) is bounded from below if φ 1 , φ 2 , φ 3 satisfies the free-leakage condition (2.2). Remark 2.2. For any ξ 1 + ξ 2 + ξ 3 = 0, there exists a constant Σ > 0 such that
if and only if the following conditions are valid:
Furthermore, the lower bound only depends on Σ 1 , Σ 2 , Σ 3 and Λ. Hence, if φ 1 , φ 2 , φ 3 satisfies the free-leakage condition (2.2), we derive ∇φ 1 + ∇φ 2 + ∇φ 3 = 0. Hence, from (2.9), we get
Remark 2.3. To form a meaningful physical system, the energy potential F(φ 1 , φ 2 , φ 3 ) defined in (2.7) must be bounded from below. For partial spreading case (Σ i > 0 , ∀i), one can assume that Λ = 0 since F 0 (φ 1 , φ 2 , φ 3 ) ≥ 0 is naturally fulfilled. For the total spreading case, Λ has to be non-zero. Moreover, in order to ensure the non-negativity of F, Λ must be large enough.
For the 3D case, it is shown in [4] that when P(φ 1 , φ 2 , φ 3 ) takes the following form, the bulk energy F is bounded from below:
where φ α (x) = 1 (1+x 2 ) α with 0 < α ≤ 8 17 . Since (2.8) is more commonly used in [2,4], we also adopt it for simplicity. It is clear that the numerical scheme developed in this paper can handle (2.9) or (2.12) without any essential difficulties.
To couple the fluid momentum into the system, the total free energy becomes
where u represents the fluid velocity.
Assuming that the fluid is incompressible and follows the generalized Fick's law, that is, the mass flux is proportional to the gradient of the chemical potential, we can derive the following three-phase Cahn-Hilliard model coupled with hydrodynamics:
)
where M is the mobility parameter, f i = ∂ i F, p is the pressure, ν is the fluid viscosity, β L is the Lagrange multiplier to ensure the free-leakage condition (2.2) and it can be derived as
The boundary conditions of the system are one of the following two types:
(i) all variables are periodic, or (ii
where n is the unit outward normal to the boundary ∂Ω . The initial conditions are given by
where the initial condition also satisfies the free-leakage condition φ 0 1 + φ 0 2 + φ 0 3 = 1.
Remark 2.4. It can be proved that the three-phase system (2.14)-(2.15) is equivalent to the following PDE system using two variables
where φ 3 and µ 3 are given by the following explicit formula:
Since the proof is quite similar to Theorem 3.1, we omit the details here.
The model equations (2.14)-(2.17) follow a dissipative energy law. By taking the L 2 inner product of (2.14) with µ i , of (2.15) with -φ it , of (2.16) with u, and performing integration by parts, we can obtain
)
Then, we take the summation of (2.24) and (2.25) for i = 1, 2, 3, and combine the obtained results with (2.26). Using (2.17) for the pressure term and (β L , (
, we obtain the energy dissipative law as
where the last inequality is derived by using (2.9) since (µ 1 , µ 2 , µ 3 ) satisfies the condition (2.23).
Remark 2.5. When deriving the PDE energy law (2.27), due to the divergence-free condition, the advection term vanishes. The advection and surface tension terms are offset by using integration by parts, i.e., ∫
This means that the advection and surface tension terms do not contribute to the total free energy or energy diffusivity, that is, they satisfy the "zero-energy-contribution" feature, which provides us with inspiration for designing a fully-decoupled scheme given in the next section.
The purpose of this section is to construct a fully-decoupled numerical scheme to solve the flow-coupled phasefield system of three components incompressible flow system (2.14)-(2.17). At the same time, considering the efficiency and accuracy of the algorithm in practice, we also expect that the scheme can satisfy linearity, secondorder time accuracy, and unconditional energy stability. The detailed process for developing such a scheme is given as follows.
First, we define the auxiliary function U (x, t) as
where B is any constant such that the radicand is always positive (Note F(φ 1 , φ 2 , φ 3 ) is always bounded from below from Remark 2.3). Second, we define another nonlocal type scalar auxiliary variable: Q(t) ≡ 1 which is written as the solution of a simple ordinary differential equation, that reads as
Then, by combining the two nonlocal variables U and Q, and the trivial evolution Eq. (3.2), the system (2.14)-(2.17) is reformulated to the following form:
where
Note (3.9) can be also written as the following form,
The variables (u, p, φ i , µ i , U, Q) constitutes a closed PDE system. The initial conditions are given as follows, Remarkably, in the new system (3.3)-(3.8), we multiply each term that satisfies the "zero-energy-contribution" characteristic with the nonlocal variable Q. It is easy to find that, after using integration by parts and the divergencefree condition, the combination of all nonlinear integrations in (3.7) is equal to zero, which means Q ≡ 1 and thus the new system (3.3)-(3.8) is equivalent to the original PDE system (2.14)-(2.17). Meanwhile, the new system (3.3)-(3.8) also holds the law of energy dissipation that can be obtained through a similar process to obtain (2.27). Since the discrete-level energy stability proof process follows the same principle, we introduce the following detailed process to make it more clear.
We multiply the L 2 inner product of (3.3) with µ i to get
We multiply the L 2 inner product of (3.4) with -φ it to get
We multiply (3.5) with 24 ϵ U to derive
By multiplying the L 2 inner product of (3.6) with u, we get
By multiplying (3.7) with Q, we get
.
(3.17)
Combining the above five Eqs. (3.13)-(3.17), using (β,
) = 0, and noting that the two terms with the same Rome numerals under curly braces cancel each other out, we derive
where the two negative terms on the right end prescribe the energy diffusive rate and
Remark 3.1. After adding the simple ODE (3.7), we can see that the derivation of the law of energy becomes slightly different from that of the original system. For example, we no longer need the term I 1 in (3.13) and II 1 in (3.16) to cancel each other. Instead, I 1 is canceled by I 2 , and II 1 is canceled by II 2 , where I 2 and II 2 are from the new ODE (3.7). In other words, for the discrete case, we can use different methods to discretize the advection (associated with I 1 ) and surface tension (associated with II 1 ), which makes it possible to design a fully-decoupled scheme.
Now, it is ready to build up a numerical scheme to discretize the new system (3.3)-(3.8) by using the second-order backward differentiation formula (BDF2). It reads as follows.
We compute ( ũ, u, p,
)
X. Yang
Computer Methods in Applied Mechanics and Engineering 376 (2021) 113589
where
S > 0 is a stabilization parameter, and the boundary conditions of the scheme are either periodic for all variables, or
In the next few remarks, we give some explanations of the scheme.
Remark 3.2. The scheme is linear and each nonlinear term is discretized using the implicit-explicit combination method. For the hydrodynamical equations, we adopt the second-order pressure-correction scheme (3.20)-(3.25)-(3.26) (cf. [44]) where ũn+1 is the intermediate velocity that follows the Dirichlet boundary conditions and the final velocity field u n+1 follows the divergence-free condition. To obtain the pressure, we just apply the divergence operator to (3.25) and then obtain the following Poisson equation for p n+1 , i.e.,
Once p n+1 is computed from (3.29), we update u n+1 by using (3.25), i.e.,
Remark 3.3. Note that the coefficient H * i actually contains explicitly processed terms f (φ 1 , φ 2 , φ 3 ), so we add the comparable stabilizer (S ∼ O(1)) in (3.22). Although this term introduces an extra error of Sδt 2 ∂ tt φ i (•), whose magnitude is comparable with the error caused by the second-order extrapolated of the nonlinear term f i (φ 1 , φ 2 , φ 3 ). In Section 4, we present enough numerical evidence to show that this stabilizer is critical to maintain the accuracy and improve the energy stability while using large time steps, see the accuracy/stability tests shown in Figs. 4.2 and 4.3.
Remark 3.4. The initialization of the second-order scheme requires all values at t = t 1 , which can be obtained by constructing the first-order scheme based on the backward Euler method. In the above second-order scheme (3.20)- (3.26), as long as we set a = 2, b = 2, c = 0, ψ * = ψ n for any variable ψ, the first-order scheme can be easily obtained.
We first prove that the discrete solutions (φ n+1 1 , φ n+1 2 , φ n+1 3 ) computed by the scheme (3.20)-(3.26) also satisfy the free-leakage condition, i.e., φ n+1
that is, there is no volume loss at all time.
Theorem 3.1. The scheme (3.21)- (3.22) is equivalent to a scheme with two variables as follows:
X. Yang
Computer Methods in Applied Mechanics and Engineering 376 (2021) 113589
with
Furthermore, from (3.34), we derive
where we use (3.33) and the definition of β * in (3.27) (that is,
. Second, we assume that Eqs. (3.21)- (3.22) are satisfied and then derive (3.31)- (3.34). We use the mathematical induction and assume that (3.33) is valid for t = t n and t = t n-1 (the validity of (3.33) at t = t 1 can be shown by performing a similar process to the first-order scheme, so the detailed proof is omitted here). For any m, we define
We take the summation of (3.21
where the advective terms vanish satisfy ∇
1 by the induction. We take the summation of (3.22) for i = 1, 2, 3 to get
Taking the L 2 inner product of (3.37) with -2δt 3 Θ n+1 , of (3.38) with C n+1 -1, and summing up the two obtained equalities, we deduce 3 4
Since all terms on the left hand side of (3.39) are non-negative, thus ∇C n+1 = 0 and ∇Θ n+1 = 0 that implies the functions C n+1 and Θ n+1 are both constants. Then (3.37) leads to C n+1 = 1 that means the (3.33) is valid. We also derive Θ n+1 = 0 from (3.38) that implies (3.34) is valid. □
Here we discuss how to implement the scheme (3.20)-(3.26) in practice. It seems that all unknown variables are nonlocally coupled together, so from the appearance, the scheme (3.20)- (3.26) does not seem to achieve the fully-decoupled structure that we expect. This reminds us that we cannot solve the scheme in any direct ways because doing so will cause a lot of time consumption. Next, we introduce the implementation process in detail, in which we make full use of the nonlocal features of the two additional variables U, Q, as shown below.
First, we use the nonlocal scalar variable Q n+1 to split (φ i , µ i , U ) n+1 into a linear combination that reads as
(3.40)
Then the scheme (3.21)-(3.22) can be rewritten as
(3.41)
According to Q n+1 , we split the system (3.41) into two sub-systems as follows,
By applying procedures similar to the second-step proof in Theorem 3.1 to the two systems (3.42) and (3.43), we derive
Note that the two subsystems (3.42) and (3.43) have the same form, so we only need to introduce the method to solve any of them, and the other follows the same line. Hence, we take the first subsystem (3.42) as an example. To solve it, we continue to use the split technique, that is, the variables (φ i,1 , µ i,1 ) n+1 are split into a linear combination form by the nonlocal variable U n+1 1 , which reads as
(3.45)
By substituting the split form of all variables in (3.45) into (3.42) and decomposing the results according to the nonlocal variable U n+1
, we obtain the following independent subsystems that read as
(3.47)
The boundary conditions for (3.46)-(3.47) are either periodic or
where
x is the explicit form. By substituting the linear form of (U, φ n+1 i ) represented by Q n+1 given in (3.40) into (3.51), we obtain
Then, according to Q n+1 , we decompose (3.52) into the following two equalities:
(3.53)
Substituting the linear form of (φ i,1 , φ i,2 ) n+1 given in (3.45) and (3.50) into (3.53), we get
(3.54)
After applying a simple factorization to (3.54), we derive
An important thing is to verify whether U n+1 1 and U n+1 2 are solvable. This can be obtained by applying simple energy estimation to the subsystem (3.47). For any ψ ∈ L 2 (Ω ) with ∫ Ω ψd x = 0, we define φ = ∆ -1 ψ to be the solution of the following Poisson equation
where the boundary condition is decided by the system (i.e., if the system is with periodic boundary condition, then ψ is periodic; if the system is with the boundary condition described in (2.19), then ∂ n φ| ∂Ω = 0.) By applying ∆ -1 to the first equation of (3.47) and combine the obtained result with the second equation of (3.47), we obtain
By taking the L 2 inner product of (3.57) with φ n+1 i,1,2 and taking the summation for i = 1, 2, 3, we derive
(3.59)
In the scheme (3.20) and (3.25)-(3.26), using (3.59) to replace the variables ( ũ, u, p) n+1 and then splitting the obtained equations according to Q n+1 , we obtain several sub-equations, which can be solved independently. More precisely, starting with (3.20), the two split variables ũn+1 1 and ũn+1 2 satisfy the following equations:
where σ 1 , σ 2 are explicit forms that are given by
We require four split variables ũn+1 i , u n+1 i , i = 1, 2 follow the boundary conditions described in (3.28), i.e., they are either periodic or satisfy:
Fourth, we solve the auxiliary variable Q n+1 . Using the split form for the variables µ n+1 i in (3.40) and ũn+1 in (3.59), one can rewrite (3.24) as the following form:
where ϑ i , i = 1, 2 are given as:
(3.65)
We need to verify (3.64) is solvable. By taking the L 2 inner product of the second equation in (3.60) with ũn+1 2 , we derive
We further take the L 2 inner product of the first equation in (3.43) with µ n+1 i,2 to get
(3.69) By using the third and fourth equality of (3.44) and applying (2.9), we derive
From the second equation of (3.53), we derive
Hence we have
Therefore, (3.66) and (3.72) imply that -ϑ 2 ≥ 0, thereby ensuring the solvability of (3.64). Finally, we update φ n+1 i , µ n+1 i for i = 1, 2, 3 and U n+1 from (3.40), ũn+1 , u n+1 , and p n+1 from (3.59). We summarize the implementation of scheme (3.20)-(3.26) as follows:
where
Proof. We multiply the inner product of (3.20) with 2δt ũn+1 in the L 2 space to get
From (3.25), for any variable v with ∇ • v = 0, we have
Hence, by multiplying 2δt 2 /3 of the above equation, we derive
By taking the inner product of (3.25) with 2δtu n+1 in the L 2 space, we have
We combine (3.75), (3.77), (3.81), and (3.82) to obtain
Computing the inner product of (3.21) with 2δtµ n+1 i in the L 2 space, we have
Computing the L 2 inner product of (3.22) with -(3φ n+1 -4φ n + φ n-1 ), we find
).
(3.85)
We multiply (3.23) with 24 ϵ U n+1 and use (3.78) to obtain 12 ϵ
(3.86)
We multiply (3.24) with 2δt Q n+1 and use (3.78) to obtain 1 2
Hence, by combining (3.83)-(3.87) and taking the summation for i = 1, 2, 3, we arrive at
where we use the following two identities
which is due to (3.33). Finally, we obtain E n+1 ≥ 0 by using (2.9) and (3.33). Likewise, we obtain (3.73) after dropping the terms in { } of (3.88) since they are positive, i.e., 3
In this section, we use the proposed algorithm (3.20)-(3.26) to perform numerical simulations, including stability/accuracy tests, 2D spinodal decomposition, 2D contact lens deposited between two stratified fluids, as well as liquid droplet rising examples in 2D and 3D. In all numerical examples, we use the rectangular computational domain. For directions with periodic boundary conditions, the Fourier-spectral method is used for discretization. For directions with boundary conditions specified in (3.28), the Legendre-Galerkin method is adopted for discretization where the inf-sup stable pair (P N , P N -2 ) is used for the velocity ( ũ and u) and pressure p, respectively, and P N is used for the phase-field variables φ 1 , φ 2 , φ 3 .
We first perform several stability tests in 2D to verify the unconditional energy stability of the fully-decoupled scheme (3.20)- (3.26). We set the 2D computational domain to be [0, 2] 2 . The initial conditions for all variables are given by The initial interface contours {φ 0 i = 1 2 } of the three phase-field variables are shown in Fig. 4.1(a). We assume that the boundary conditions are periodic and space is discretized by using the Fourier-spectral method where 257 Fourier modes are used to discretize each direction. Using the above initial conditions and a time step of δt = 0.01 2 3 , we implement the algorithm until a steady-state solution is obtained, and plot the fluid interface with the velocity field at the steady-state in Fig. 4.1(b). We can see that due to the balanced surface tension coefficients of the three fluid components and coarsening effects, the two droplets squeeze together to form a contact angle of 120 • . In Fig. 4.1(c), we plot the profile of the pressure p at the steady-state using a color diagram.
Next, we verify whether the scheme maintains energy stability unconditionally using any time step. For convenience, we represent the developed decoupled scheme (3.20)-(3.26) using scalar auxiliary variables by DSAV. Meanwhile, in order to illustrate the advantages of the developed decoupling technique and the added stabilizer in improving energy stability, we also test the stability of the scheme DSAV, but the variable Q n+1 and S are removed (i.e., assume Q n+1 ≡ 1 and S = 0). For convenience, we refer to this version as SAV. In Fig. 4.2(a), we plot the evolution curves of the total free energy (3.74) calculated by the scheme DSAV using different time steps. We find that all calculated energy curves show monotonic attenuation, which confirms the unconditional stability of DSAV. When the time step is relatively large, the energy curve with a large time step has an observable deviation from the energy curve with a small time step. This is because the errors obtained using large time steps are also large. When the time step is less than 0.01 2 3 , the five energy curves obtained are very close, which is why we use the time step size 0.01 2 3 to obtain the equilibrium solution shown in Fig. 4.1. For comparison, in Fig. 4.2(b), we plot the energy evolution curves calculated by SAV, and observe that it blows up while a large time step size is used, and only decays when δt ≤ 0.01 2 5 . Through performing mesh refinement tests in time for the example above, we further test the convergence order of the developed scheme DSAV. Since the exact solutions are unknown, the numerical solution obtained by the scheme DSAV using a very small time step (δt = 1e-9) is regarded as an exact solution for calculating the approximate error. Then, we plot the L 2 errors of all variables at t = 0.4 obtained by changing the time step size from 0.01 to 0.01 2 7 with a factor of 1/2. The convergence rate is shown in Fig. 4.3(a), where we find that the scheme DSAV always exhibits almost perfect second-order accuracy. In Fig. 4.3(b), we compare the accuracy of DSAV and SAV by plotting the arithmetic mean of L 2 numerical errors of the three phase-field variables. When the time step size is large (δt > 0.01 2 6 ), SAV completely loses the convergence order. When the time step size is small (δt ≤ 0.01 2 6 ), it has the second-order accuracy. For comparison, DSAV consistently shows good second-order accuracy in all tested time steps. If we compare the size of the error when δt ≤ 0.01 2 6 , we find that the error obtained by SAV is smaller than that obtained by DSAV. This is because the stabilization terms increase some splitting errors.
In this example, we study the equilibrium pattern obtained after phase separation (or spinodal decomposition). The initial conditions are set as a homogeneous ternary mixture but with some small perturbations, which read as Therefore, we adjust the three surface tension parameters (σ 12 , σ 13 , σ 23 ) to investigate whether the contact angles under steady-state follows the theoretical predictions from (4.8). We use four partial spreading cases and two total spreading cases in the computations. Theoretical prediction values of the contact angles according to the given surface tension parameters are shown in Table 4.1. In Fig. 4.7, We first set the three surface tension parameters to (σ 12 , σ 13 , σ 23 ) = 10(1, 1, 1) and use two different gravity parameters g 0 = 90 and 100. Figs. 4.10 and 4.11 shown numerical solutions at various times, of which snapshots of the profiles of 1 2 φ 1 + φ 2 are plotted. When the gravity constant is relatively weak (g 0 = 90), we find that the droplet is captured by the interface of the two layered fluids. When the gravity constant is large (g 0 = 100), the droplet passes through the interface and enters the upper fluid.
Furthermore, we change the surface tension parameter to (σ 12 , σ 13 , σ 23 ) = 10(1, 0.6, 0.6), and also use two different gravity parameters g 0 = 90 and 100, shown in Figs. 4.12 and 4.13. Similar to the previous two simulations, similar phenomena are observed, including the low gravity parameter leading to capture and high gravity parameter leading to interface penetration. In addition, when g 0 = 90, we see that the contact angle caused by the suspension In Fig. 4.14(a) and (b), we adopt the surface tension parameter (σ 12 , σ 13 , σ 23 ) = 2(1, 0.6, 0.6) and different gravity parameters g 0 = 40 and g 0 = 50, respectively. We use different colors to plot the isosurfaces of {φ 1 = 0.5} (yellow) and {φ 3 = 0.5} (red). Similar to the 2D simulation, lower gravity causes the droplet to be captured by the interface, while higher gravity cause the droplet to penetrate the interface. In addition to the capture/penetration phenomenon, we also find that when the gravity constant is large, a long filament forms after the droplet penetrates the interface, and then the rupture of the filament occurs (at t = 5.4).
In this paper, with the help of two nonlocal auxiliary variables, a novel second-order fully-decoupled numerical algorithm is developed to solve the highly nonlinear phase-field model of three-phase incompressible flow. The scheme only needs to solve several decoupled linear elliptic equations with constant coefficients at each time step to obtain a numerical solution with second-order time accuracy. Solvability and unconditional energy stability have also been rigorously proven, and a large number of numerical simulations have been performed in 2D and 3D to show the accuracy and stability of the scheme numerically. Moreover, the novel decoupling method proposed in this paper can be widely applied to a large number of coupling models, as long as the nonlinear coupling terms satisfy the so-called "zero-energy-contribution" characteristic.