Highly conservative Allen–Cahn-type multi-phase-field model and evaluation of its accuracy

In the engineering field, it is necessary to construct a numerical model that can reproduce multiphase flows containing three or more phases with high accuracy. In our previous study, by extending the conservative Allen–Cahn (CAC) model, which is computationally considerably more efficient than the conventional Cahn–Hilliard (CH) model, to the multiphase flow problem with three or more phases, we developed the conservative Allen–Cahn type multi-phase-field (CAC–MPF) model. In this study, we newly construct the improved CAC–MPF model by modifying the Lagrange multiplier term of the previous CAC–MPF model to a conservative form. The accuracy of the improved CAC–MPF model is evaluated through a comparison of five models: three CAC–MPF models and two CH–MPF models. The results indicate that the improved CAC–MPF model can accurately and efficiently perform simulations of multiphase flows with three or more phases while maintaining the same level of volume conservation as the CH model. We expect that the improved CAC–MPF model will be applied to various engineering problems with multiphase flows with high accuracy.

Eulerian numerical methods are typically used for expressing multiphase flows; these include the volumeof-fluid (VOF) method [24,40,45], the level-set (LS) method [11,38,41], and the phase-field (PF) method [2,12,13,19,[33][34][35].In the VOF method, the flow is represented by assigning the volume fraction of a phase to each numerical grid according to the Eulerian method, which has excellent volume conservation.However, the interface between two phases becomes discontinuous, and a complex numerical algorithm must be solved to make the interface continuous [24].In the LS method, the interface is linearly expressed by a distance function called the LS function.When solving this function, it is necessary to reinitialize it each time step so as to satisfy the governing equation of the fluid and the law of volume conservation, which incurs a high calculation cost; moreover, the volume error increases with time [28,41].In the PF method, the interface is treated as a smoothly varying region of PF variables.No special treatment is required to maintain and capture the interface, and the evolution of the interface can be represented by solving the time-evolving equation for the PF variables defined as continuous values.Furthermore, the law of volume conservation, which is important in multiphase flow simulations, is strictly satisfied by solving the Cahn-Hilliard (CH) equation [9].The high computational cost is a drawback, but if this problem is solved, the PF method can be considered the most promising model to represent multiphase flows.Multiphase flows of three or more phases can easily be handled by increasing the number of PF variables accordingly; this is known as the multi-phase-field (MPF) method [1,2,5,12,13,16,19].
The CH-MPF model [5,16,19], which is an extension of the CH equation, has been primarily developed as an MPF model for multiphase flows.Although the CH equation satisfies the law of volume conservation because of the inclusion of the equation of continuity, it is a fourth-order parabolic partial differential equation and requires a very small time increment, which results in a high computational cost.To address this problem, Beckermann et al. [31] considered subtracting the curvature effect from the traditional Allen-Cahn (AC) equation [8], which is non-conservative.However, this model is not strictly conservative; thus, the conservative Allen-Cahn (CAC) equation was later proposed by Chiu et al. [10].The CAC equation is a second-order differential equation, and simulations using the CAC model can be run stably with a significantly reduced computational cost because of the larger time increment in comparison to the CH equation in the case of using an explicit method.Furthermore, it has been indicated that the CH equation includes the curvature effect of the interface; therefore, when coupled with the Navier-Stokes equation considering surface tension, the surface tension is considered doubly, and this effect cannot be properly represented [34].Thus, in multiphase flows with three or more phases, CAC-MPF models have been actively studied [1,2,12,13,43].
Hu et al. [12] and Abadi et al. [1] coupled the CAC-MPF equation with the lattice Boltzmann equation.Huang et al. constructed a CAC-MPF model with bounded mapping, enabling two-dimensional three-phase fluid calculations while satisfying volume conservation in each phase [13].Yang and He [43] developed a finite element scheme to solve a CAC-MPF model for a system consisting of three phases while maintaining unconditional energy stability.Moreover, Aihara et al. constructed a CAC-MPF model for multiphase flows with an arbitrary number of phases using the Lagrange multiplier [2].The CAC-MPF model is expected to be applied to three-dimensional simulations and more complicated multiphase flow phenomena because of its lower computational cost and better morphology-preserving property, when compared to the CH-MPF model.However, the conservativeness of the CAC-MPF model strongly relies on the modeling method.A general CAC-MPF model is constructed by adding a Lagrange multiplier [26] to the CAC equation representing two-phase flow, and its conservation is determined by the formulation of the multiplier.
In this study, we newly construct a modified CAC-MPF model by extending our previously constructed CAC-MPF model [2] to a more conservative system.Then, the validity of the new model is evaluated through a comparison with five models: three CAC-MPF models and two CH-MPF models.The rest of this paper is organized as follows.Section 2 derives the proposed CAC-MPF model and also show other CAC-MPF and CH-MPF models for comparison.Section 3 presents the results of five kinds of two-dimensional simulations to demonstrate the effectiveness of the proposed model.Other MPF models for multiphase flows are also discussed.Finally, the conclusions are presented in Sect. 4.

Multi-phase-field equations
In this section, we derive the CAC-MPF model proposed in this study.We also describe the CH-MPF model proposed by Kim et al. [19] and the model from which the curvature effect is removed.
The free energy F considering only the interfacial energy in a two-phase system with a single PF variable φ can be expressed, for example, as follows [36].
Note that φ =1 and 0 represent the PF variable inside both phases, and φ is defined to change smoothly from 0 to 1 at the interface.The first term in the right-hand side of Eq. ( 1) denotes the gradient energy density, and the second term denotes the double well potential, where ais the gradient coefficient and W is the height of the energy barrier.Substituting Eq. ( 1) into the AC equation ∂φ /∂t = -M AC δ F/δφ yields the following AC equation [8]: where M AC is the PF mobility in the AC equation.Considering the one-dimensional equilibrium state in Eq.
(2), we have In Eq. ( 3), the thickness of the interface is defined as a region where λ < φ < 1-λ, which determines the constant b as b = 2tanh −1 (1-2λ).Setting the interfacial energy σ and the thickness of the interface δ, the gradient coefficient a and energy barrier height W are expressed by the following relations [36]: The CAC equation, which is the conservative form of Eq. ( 2), is expressed as follows [10]: Next, the MPF equation is derived by extending Eq. (5) to a system with an arbitrary number of phases, where multiple PF variables φ i are introduced to represent multiple phases in the MPF model.The variable φ i takes 1 in phase i and 0 in the other phases and smoothly changes from 0 to 1 in the interfacial region.The PF variables are not independent and must satisfy the following constraints: Here, N is the number of phases in the entire domain.Differentiating both sides of Eq. ( 6) with time, we obtain The free energy functional F is established using the multiple PF variables φ i defined earlier.By simply superposing N equations of Eq. ( 1) and expressing F, the following equation is obtained: where a and W are assumed to be constant in all interfaces.To satisfy Eq. ( 6), the Lagrange multiplier is applied to express the free energy functional Fof the entire system: Here, α is the Lagrange multiplier.The functional derivative of Eq. ( 9) with φ i is Following the same procedure as the derivation of Eq. ( 5), we obtain the following CAC-MPF equation: Note that Γ AC = M AC a 2 is the mobility corresponding to the diffusion coefficient.Additionally, the constant B is B = √ 2W a = 2b/δ, which corresponds to the inverse of the thickness of the interface.The last term in Eq. (11) indicates the effect of multiphase interactions.The Lagrange multipliers must be determined in such a way that each PF variable can exist stably at multiple points while satisfying Eq. (6).
Substituting the functional in Eq. ( 10) into the AC equation ∂φ i /∂t = -M AC δ F/δφ i and further into the constraint Eq. ( 7), we obtain We then apply the Lagrange multiplier [16] proposed by Lee et al. in the following equation: As shown here, in the formulation by Lee et al., different multipliers are adopted for each phase, where qis a constant.Lee et al. have shown that as the value of q increases, the influence of energy stabilization among multiple phases because of the inflow of additional phases on multiple points can be suppressed [16].In this study, q = 2 was used based on previous results [2,16,19].In the functional derivative of Eq. ( 13), the first term a 2 ∇ 2 φ i on the right-hand side of Eq. ( 10) is omitted because it is zero when summed.By assuming α = α i and substituting Eq. (13) into Eq.( 11), we can derive Equation (14) expresses that the summation of the differential forms of the double well potential term in Eq. (1), which prevents instability because of energy dissipation at higher-order points, is simply added to the CAC equation to express a two-phase system as a weight function.
Substituting Eq. ( 11) into Eq.( 7) yields the following equation for the Lagrange multiplier: As in the derivation of Eq. ( 13), we follow the formulation of Lee et al. [19] and use the following multipliers, which are different for each phase: Unlike Eq. ( 14), the multiplier in Eq. ( 16) includes the spatial derivative to represent the interaction between phases while considering the diffusion of the phases at the higher-order points.In the same way as in Eq. ( 14), by substituting Eq. ( 16) into Eq.( 11) as α = α i , we obtain Eq. ( 17) is the model constructed in our previous study [2].Although our previous study confirmed that Eq. ( 17) had sufficient accuracy, the Lagrange multiplier shown in Eq. ( 16) was not in a strictly conservative form, and there is a possibility that the phase volume is not conserved at a higher-order point.Therefore, Eq. ( 16) should be modified to a more accurate form.For this reason, Eq. ( 16) is modified such that the sum of the right-hand side becomes a conservative formulation while satisfying Eq. ( 12), and we obtain By substituting Eq. ( 18) into Eq.( 11), we can derive the improved CAC-MPF model developed in this study: The MPF equations shown in Eqs. ( 14), (17), and ( 19) can express the time evolution of each phase while satisfying Eq. ( 6) in the MPF model.
Next, the MPF equations based on the CH equation, or CH-MPF equations, [9] are derived in the same way as the CAC-MPF equations, which are used in a comparison with the CAC-MPF equation in the next section.Substituting the free energy F in a two-phase system with a single PF variable shown in Eq. ( 1) into the CH equation ∂φ /∂t = ∇[M C H (∇δ F/δφ]], we obtain the following time evolution equation: where Γ C H = M C H a 2 is the mobility in the CH equation with the PF mobility M C H .To extend Eq. ( 20) to the MPF equation, the same procedure used to derive Eq. ( 14) from Eq. ( 6) is performed for Eq. ( 20), yielding [19] As the CH-MPF equation shown in Eq. ( 21) includes the curvature effect, the new CH-MPF model is derived by removing the curvature from the diffusion term as follows: In the next section, to demonstrate the usefulness of Eq. ( 19) derived in this study, we compare it with Eqs. ( 14), ( 17), (21), and (22).For ease of notation, we denote the models obtained using each equation as follows: Eq. ( 14): CAC1 model; Eq. ( 17): CAC2 model; Eq. ( 19): CAC3 model; Eq. ( 21): CH1 model; and Eq. ( 22): CH2 model.Note that when these equations are used under advection conditions, the advection term ∇(φ i u ) is added to each equation.
The CAC3 model, Eq. ( 19), is similar to the N -phase Conservative Level Set Free-Energy-based Surface Tension model (NCLS-FEST) by Howard and Tartakovsky [11].The difference is with and without the term with Δφ j in the summation, which was caused by the base equations from which the derivation was started.The CAC3 model was derived starting from the free energy functional, Eq. ( 1), and the Allen-Cahn equation, while the NCLS-FEST was derived based on the conservative level set method [22].

Governing equations for flow
Assuming an incompressible flow, the flow velocity u is obtained by solving the Navier-Stokes equations, taking into account the equation of continuity and the surface tension force at the interface between the different phases i and j, as shown in the following equation: where ρ is the density, p is the pressure, μ is the viscosity, g is the gravitational acceleration, and e g is a unit vector indicating the direction of gravity.The density and viscosity coefficients are expressed as ρ = Σ i ρ i φ i and μ = Σ i μ i φ i , respectively, where ρ i and μ i are the density and viscosity in each phase expressed by the PF variable φ i .In addition, SF indicates a volume force vector because of the surface tension and is expressed as by adding up the surface tension vectors for each interface.Here, γ i j is the interfacial energy at the i-j interface; sf i and sf j are the curvature vectors due to the gradients of φ i and φ j , and ω(φ i , φ j ) is a function that determines the contribution range of the surface tension in the interfacial area.The curvature vector sf i has been modeled by Blackbill [7] for a two-phase flow using which is expressed by the inner product of the curvature κ and the unit vector n in the normal direction.[16], and the integral over the entire region satisfies which results in the constant E = −15δ /b.

Nondimensionalization
The MPF equations and governing equations of the fluid are nondimensionalized using the following relations: where L is the representative length, U is the representative velocity, ρ re f is the representative density, μ re f is the representative viscosity coefficient, Pe AC and Pe C H are the Peclet numbers, and the letters with "-" above them represent dimensionless quantities.Using Eq. ( 28) to nondimensionalize the CAC1, CAC2, CAC3, CH1, and CH2 models with the added advection term, we obtain the following equations.CAC1 model: CAC2 model: CAC3 model: CH1 model: CH2 model: As Pe AC and Pe C H are parameters to recover the phase-field profile to the equilibrium in one step, they rely on the maximum flow velocity at the interface, | ūmax | [10].The nondimensionalization forms of Eqs. ( 23) and ( 24) are expressed as follows: where the Reynolds number Re, Weber number We i j , and Froude number Fr at the i − j interface are, respectively, expressed as follows:
The simplified marker-and-cell (SMAC) method [20] is used to solve the Navier-Stokes equations assuming an incompressible fluid.The pressure increment is computed by solving the Poisson equation derived from the equation of continuity (Eq.( 34)) and the two-dimensional discretized Navier-Stokes equation (Eq.( 35)) using the intermediate velocity of n and n+ 1 steps.The Poisson equation is solved using the geometric multigrid (GMG) method, and the velocity u n+1 and pressure p n+1 that satisfy the equation of continuity in n + 1 steps are obtained from the pressure increment.In this study, the advection, diffusion, and surface tension terms in the Navier-Stokes equations are discretized using the fifth-order WENO scheme, second-order central difference method, and fourth-order accurate central difference method, respectively.The efficiency of the calculation is improved by introducing the preconditioned GMG-BiCGStab [39].We developed the program code using CUDA C, and the simulations were carried out using an NVIDIA Tesla P100 GPU.
The solving procedure of Eqs. ( 29)-( 35) is as follows: (1) Set the initial condition and initialize all variables φ 0 i , ū0 , and p0 .Here, a superscript indicates a time step.(2) The PF variable φ n+1 i 0 is calculated by solving each MPF model, where a superscript of the curly bracket represents a stage of the RK method and R H S n i corresponds to the right-hand side of Eqs. ( 29)- (33).In the calculation of R H S n i , the second-order central difference method is uniformly used for the Laplacian calculation.
In the NCLS-FEST [11], the equation corresponding to Eqs. ( 29)-( 31) was solved to steady state.This can recover the variable to an equilibrium profile with high accuracy.On the other hand, the computational cost becomes larger, and, thus, the explicit scheme employed in the present study is more efficient and versatile.

Simulations
The accuracy of the developed CAC-MPF model is evaluated via comparison of five MPF models in three types of two-dimensional multiphase flow simulations.First, the Ostwald ripening is investigated to assess interfacial morphology preservation for each model.Next, advection simulations of rotating circular phases are performed to evaluate volume conservation.Then, pressure jump and liquid lens tests are conducted, and the results are compared with theoretical values and previous studies.Finally, we simulate the bubble rising through the liquid-liquid interfaces.

is expressed by
where r is the coordinate to the direction of interface normal with r = 0 at φ eq i = 1/2.Equation ( 44 29)- (33) without the advection term.To observe the morphological changes of the phases in detail, the simulations were performed until t = 1835, where there were no more changes in the interfacial morphology.
Figure 2 indicates the distribution of φ i in a log notation on the A -A and B -B lines (which are indicated by dashed lines in Fig. 1).As the three CAC and CH2 models gave exactly the same results, only the results for CAC3 are shown in Fig. 2a. Figure 2b shows the results for the CH1 model.The dashed line shows the initial distribution, and the solid lines show the distribution at t = 1835, where blue, gray, and black lines indicate phases 1, 2, and 3, respectively.As shown in Fig. 2a, the solid and dashed lines around φ i > 0.01 almost overlap, indicating that the four circular phases retain their initial morphology.Conversely, the results in Fig. 2b reveal that the small circular phases 1 and 2 are absorbed by the larger ones, resulting in Ostwald ripening, in which the larger circular phases become even larger.This is because the CAC and CH2 models exclude the curvature effect, whereas the CH1 model includes the curvature drive [18,34].These results demonstrate that it is important to remove the curvature drive from the PF equation when using the PF method as an interface tracking method [34].Figure 2(a) also indicates that the PF distribution changes in the range φ i ≤ 10 −2 , because the PF variable has a hyperbolic tangent distribution as shown in Eq. (39).Although the distribution of each phase changes slightly, the interfacial morphology is retained, implying that there is no issue with morphological retention.In addition, in this evaluation, the area error with respect to the initial state is within the order of 10 −12 for all models, indicating that there is no issue with volume conservation.

Rotating disc with two phases
We evaluate the volume conservation of the MPF models considering advection.As shown in Fig. 3, we consider, within a 5L × 5L (L = 50Δx) computational domain, a circular phase 1 with radius L/2 (φ 1 = 1, φ 2 = φ 3 = 0), a crescent-shaped phase 2 (φ 2 = 1, φ 3 = φ 1 = 0) bordering it, and remaining phase 3 (φ 3 = 1, φ 1 = φ 2 = 0).The initial positions of the circle centers in phases 1 and 2 are (5L/2, 15L/4) and (5L/2, 13L/4), respectively.For this state, a constant angular velocity is applied so that phases 1 and 2 rotate clockwise, as shown by the dotted line in Fig. 3, with the center of rotation being at the center of the domain (5L/2, 5L/2).Note that the time-evolving equations of Eqs. ( 29)- (33) for φ i are solved to evaluate the effect of the advection term.The coefficients in the time evolution equations were determined as B = 1, Pe AC = 8.85× 10 2 , and Pe C H = 8.85× 10 3 .The time increments for the CAC and CH models were set as Δt AC = Pe AC Δ x2 /4 = 2.51× 10 −3 and Δt C H = Pe C H Δ x4 /50 = 7.84× 10 −6 , respectively.The time increment for the CAC model was 320 times larger than that for the CH model.Thus, 2,500Δt AC for the CAC model and 800,000Δt C H for the CH model were required for phases 1 and 2 to rotate one cycle, and simulations of 12 cycles were performed.Here, the φ i is set to the zero-Neumann condition at all boundaries.
The interfacial contour lines are shown in Fig. 4 as Σφ 2 i = 0.58, 0.68, and 0.82 after 12 cycles, where the colored lines are the results after 12 cycles and the initial state is indicated by the solid black line.The differences between the models cannot be confirmed in the scale of Fig. 4. Thus, the results for the CAC3 model are displayed as representative.From Fig. 4, it can be observed that the interface thickness at 12 cycles becomes larger than the initial one.Figure 5 shows a close-up view of the square area indicated by the solid black line in Fig. 4. Figure 6 represents the time variation of the area in phases 1 and 2 with respect to the initial value, in a form of Error = (A i -A ini )/A ini .The area of phase i is computed by A i = ∫ Ω φ i dΩ, where Ω is the computational domain.In Fig. 4, the results for the CAC1, CAC2, CAC3, CH1, and CH2 models are shown by solid red, blue, peach, yellow, and green lines, respectively, with circles and triangles representing the volume error for areas A 1 and A 2 for φ 1 and φ 2 , respectively.Figure 4 shows that the interface thickness is slightly greater than the initial one because of the advection term; however, the morphology of each phase is stable even after 12 cycles.At the triple point, there is considerable deviation from the initial interfacial morphology, and the degree of this deviation is different for each model.Figure 5a-c indicate that the deviation from the initial after 12 cycles is larger in the order of CAC3, CAC2, and CAC1, and the accuracy is higher in the order of CAC1, CAC2, and CAC3.By contrast, according to Fig. 5d and e, the morphology for the CH1 and CH2 models is almost consistent with that of the CAC3 model.Figure 6 shows that the error of the CAC1 and CAC2 models is in the order of 10 −8 , while that of the CAC3 and two CH models is in the order of 10 −14 , indicating that volume conservation of the CAC3 is highly accurate in the same order with the CH models.Figure 6 also confirmed the fact that the errors of A 1 and A 2 are almost identical for the CAC3 and CH models, whereas the error of A 2 with a large curvature is larger for the CAC1 and CAC2 models.Furthermore, in Fig. 6, the fluctuations in the curve are smaller for the CAC3 model than for the CH models.
Therefore, the CAC3 model developed in this study has the same level of volume conservation as the CH model, and the simulation is stable.Moreover, the time increments are much larger than those of the CH model, which makes it possible to perform multiphase flow simulations with high accuracy and speed.

Pressure jumps in static droplets
The pressure jumps are evaluated in a system with statically placed three bubbles with different surface tensions [11,16].Figure 7 shows the calculation domain and initial conditions for the pressure jump test.In a computational region of 11L × 11L (L = 64Δx ), three circular phases with radius L are placed.The blue phase 1 (φ 1 = 1, φ 2 = φ 3 = φ 4 = 0), the gray phase 2 (φ 2 = 1, φ 3 = φ 4 = φ 1 = 0), and the yellow phase 3 (φ 3 = 1, φ 4 = φ 1 = φ 2 = 0) are centered at point C 1 (33L/4, 33L/4), C 2 (11L/2, 11L/2), and C 3 (11L/4, 11L/4) respectively, and the other phase arranged to the white phase 4 (φ 4 = 1, φ 1 = φ 2 = φ 3 = 0).We set the non-gravitational acceleration vector g = (0, 0), representative length L = 1 m, representative velocity U = 1 m/s, representative density ρ re f = 10 kg/m 3 , representative viscosity μ re f = 1 Pas, and dimensionless constants for each phase as follows: ρ1 = ρ2 = ρ3 = ρ4 = 1, μ1 = μ2 = μ3 = μ4 = 1, B = 70.4,Pe AC = Pe C H = 1.61,Re = 10, We 12 = We 13 = We 23 = 10, We 14 = 2, We 24 = 4, and We 34 = 1.The Weber number (We 14 , We 24 , We 34 ) = (2, 4, 1) corresponds to the interfacial energies (σ 14 , σ 24 , σ 34 ) = (5 J/m 2 , 2.5 J/m 2 , 10 J/m 2 ).The Peclet numbers were determined as the interfacial morphologies and are appropriately kept in the liquid flow.The time increments for the CAC and CH models were determined as Δt AC = 9.81× 10 −5 and Δt C H = 1.92× 10 −9 .Here, all boundary conditions for the PF variable φ i , the velocity u, and the pressure Figure 8 indicates the pressure profile along the A -A line in Fig. 7 for the CAC3 model.Because the differences between five models are not visible in Fig. 8, only the results of the CAC3 model are shown as representative.The dashed lines in Fig. 8 are the theoretical values of pressure in each phase, P theor y = σ i j /L.From Fig. 8, it can be observed that the pressures in each circular phase are in good agreement with the theoretical values.The sharp drop of the pressure within the interface is because of a numerical fluctuation of a Fig. 7 Initial settings and computational conditions for circular phases 1, 2, and 3 placed in the phase 4 Fig. 8 Pressure profiles along the line A − A in Fig. 7.The dashed lines are the Laplace pressure in each phase gradient in the distribution, and the entire trend of pressure distributions matches with that in previous studies [11].Furthermore, the accuracy of pressure inside each phase, P calc , for the theoretical value is estimated as an error defined by Error = (P calc − P theor y )/ P theor .The P calc is calculated as the averaged value of grid points satisfying φ i > 0.999.The calculated errors are 1.21× 10 −3 , 1.30× 10 −3 , and 1.13× 10 −3 for phase 1, 2, and 3, respectively, and these were almost the same among the five MPF models.It can be confirmed that the errors are sufficiently small.These results confirm that all the MPF models can accurately represent pressure jump among multiple phases, and the results are in good agreement with other models [11,16].

Liquid lens test
We evaluate an equilibrium morphology of liquid lens placed on a liquid-liquid interface.Figure 9 shows the computational region 5L × 5L (L = 128Δx ) and initial profiles of three phases.The blue circular phase 1 (φ 1 = 1, ϕ 2 = φ 3 = 0) with radius L is located at (5L/2, 5L/2), and the white phase 2 (φ 2 = 1, φ 3 = φ 1 = 0) is placed as the upper half of the domain separated by a central line y = 5L/2.The other gray phase 3 (φ 3 = 1, φ 1 = φ 2 = Figure 10 shows the interfacial morphologies at the steady state calculated using the CAC3, where the lines indicate a region with Σφ 2 i < 0.82.As shown in Fig. 10, the lens morphology becomes thinner as the Weber number on the 1-2 and 2-3 interfaces increases, that is, as the surface tensions decrease.Table 1 lists the lens width d calc of the phase 1, which corresponds to d in Fig. 9, computed for the five MPF models, where the values are shown as the errors against the theoretical width d theor y [30] The results of Ref. [11] are also included in Table 1.Here, the value d calc is obtained as the distance between two points satisfying φ 1 > 0.3, φ 2 > 0.3, and φ 3 > 0.3.From Table 1, it is confirmed that the d calc agree well with the d theor y with high precision within an error of 1× 10 −2 order.Moreover, it can be shown that the results using the three CAC models are generally more accurate in the order of CAC1, CAC2, and CAC3, and the errors between CAC3, CH1, and CH2 are almost the same.Compared with the results of [11], the accuracy of the CAC3, CH1, and CH2 models is better.
From the above results, it is shown that CAC3, CH1, and CH2 models are equally accurate and enable to express lens morphology with a high accuracy.

Bubble rising in three phases
Simulations of a bubble rising in three phases with different densities [2,19] are performed to verify the performance of each MPF model for the problem with liquid flow.The computational domain L × 4L (L = 64Δx) and initial conditions are shown in Fig. 11; phase 1 (φ    i < 0.82.The morphological changes of the bubble were similar to Fig. 12 for other MPF models.From Fig. 12, we can clearly see that the bubble morphologies floating up through balk liquid and the liquid-liquid interface differ depending on the position because of the differences in the buoyancy force and interface energy.
Figure 13 shows the temporal changes of volume error of the area of a rising bubble, A 1 , with respect to the initial volume.Similar to the results shown in Fig. 6, the errors of the CAC1 and CAC2 models are larger than those of the other models.A comparison of the CH1 and CH2 models indicates that both errors are small; however, the CH1 model has an error that is one-order larger compared to that of the CH2 model.In addition, the CAC3 model initially shows a one-order larger error than the CH2 model, whereas this error is almost constant during the simulation.On the other hand, the error of the CH2 model increases during the simulation, and finally, the error becomes larger than that of the CAC3 model.In summary, the volume conservation at the end of the simulation was higher for the CAC1, CAC2, CH1, CH2, and CAC3 models, in that order.Thus, it can be seen that the CAC3 model can satisfy the volume constraint accurately and stably.
Figure 14 shows the velocity of the rising bubble depending on the y-position of the bubble's center of gravity.As indicated in Fig. 14, the rising velocity of the bubble is the smallest at two liquid-liquid interface positions, whereas those are slightly higher than the initial position of the two interfaces at y = 4L/5 and 5L/2.At y = 4/5, the bubble velocities of all models are almost the same, whereas at y = 1, the bubble velocities are larger for CAC1, CAC2, and the other models in that order.The difference in velocity at the same y-coordinate increases as the bubble rises, particularly after passing through the liquid-liquid interface.The results of the CAC3 and two CH models are found to be in close agreement throughout the simulations.
These results confirm that the CAC3 model shows stable and accurate volume conservation in the evaluation of bubble dynamics with liquid phase flows.

Conclusions
In this study, we newly constructed an improved CAC-MPF model, in which the Lagrange multiplier term of the CAC-MPF model constructed in our previous study was improved to a conservative system.The accuracy of In the future, we plan to deal with actual multiphase flow problems with more than three phases by extending the CAC-MPF model to a three-dimensional problem.As in the PF method, the time evolution equations are same independent on the dimension, the CAC-MPF equation, Eq. ( 31), can be directly used in the theredimensional simulations.Nevertheless, this model has a problem in that the computational cost increases in proportion to the number of phases because all the phases to be considered have some values in the entire domain.In the future, we will address this problem and construct a model that can treat an arbitrary number of phases in a short time.

( 6 ) Set φ n+1 i 0 ,, 3 ,
ūn+1 0 , and pn+1 0 in steps(1) to(5) as the initial stage of the time derivative values.By solving the RK method, φ n+1 i k ūn+1 k , and pn+1 k are approximated precisely by using the values of the previous stage, where k is the number of a stage of approximation.The next value of the stage k+ 1 is determined by the present one of the stage k.In the third-order RK method, the final values φ n+1 i ūn+1 3 , and pn+1 3 are the approximations and finally updated as the next time ones.(7)The step number is replaced as n = n + 1.
) is dimensionless and the MPF version of Eq. (3).The dimensionless constants in the MPF equations are set to B = 1 and Pe AC = Pe C H = 5.46.The time increments are set to Δt AC = Pe AC Δ x2 /4 = 4.58× 10 −2 and Δt C H = Pe C H Δ x4 /50 = 3.67× 10 −3 .The zero-Neumann condition for φ i is set at all boundaries.The morphological change from Fig. 1 is simulated by solving the MPF equations in Eqs. (

Fig. 2 Fig. 3 Fig. 4
Fig. 2 Profiles of φ 1 (blue line), φ 2 (gray line), and φ 3 (black line) along the lines of A -A and B -B in Fig. 1 at the initial (dashed line) and final step at t = 18350 (solid line) (color figure online)

Fig.Fig. 6
Fig. Close-up views around the triple point enclosed by the square in Fig. 4. In addition to the lines shown in Fig. 4, the results at t = 24π for the CAC3 model (red line) are drawn for all figures except Fig. 5(c) (color figure online)

Fig. 9
Fig. 9 settings and computational conditions in the liquid lens tests

Fig. 11
Fig.11settings and computational conditions for the bubble (phase 1) rising through three liquid phases and two liquidliquid interfaces

Fig. 12 1 Fig. 14
Fig. 12 Snapshots of temporal bubble morphological changes in case of the CAC3 model.The nondimensional time of each snap is t = 0, 2.34, 5.46, 7.03, 10.15, 13.28, 16.40, 17.18, 18.75, and 23.43 from left to right.The color indicates the phases in the same way as in Fig. 11, and the solid black lines represent the liquid-liquid and gas-liquid interfaces

Table 1
Accuracy of lens width d (45) .Values are expressed as error defined as (d calcd theor y )/d theor y with the theoretical value d theor y by Eq.(45)