Mass-conservation-improved phase field methods for turbulent multiphase flow simulation

The phase field method has emerged as a powerful tool for the simulation of multiphase flow. The method has great potential for further developments and applications: it has a sound physical basis, and when associated with a highly refined grid, physics is accurately rendered. However, in many cases, especially when dealing with turbulent flows, the available computational resources do not allow for a complete resolution of the interfacial phenomena and some undesired effects such as shrinkage, coarsening and misrepresentation of surface tension forces and thermo-physical properties can affect the accuracy of the simulations. In this paper, we present two improved phase field method formulations (profile-corrected and flux-corrected), specifically developed to overcome the previously mentioned drawbacks, and we benchmark their performance versus the classic one. The formulations are first tested considering the rise of a bubble in a quiescent fluid and the interaction of two droplets in laminar shear flow; then, their performances are compared in the simulation of a droplet-laden turbulent flow. The aim of this work is to review and benchmark the different phase field method formulations, with the final goal of laying down useful guidelines for the accurate simulation of turbulent multiphase flow with the phase field method.


Introduction
The accurate simulation of turbulent multiphase flow is of crucial importance in a wide range of applications, from raindrop formation [11] to the pharmaceutical industry [15]. The description of such flows poses several challenges, among them the description of topological changes, viscosity/density contrasts and surface tension forces [17]. In this context, the phase field method emerges as a powerful alternative to the more conventional methods such as front tracking (FT) [41], volume of fluid (VOF) [18,38] and level set (LS) [29,30]. In recent years, the phase field method has been used to simulate turbulent multiphase flow [1,2,6,[33][34][35][36][37], binary alloys [12,42,43] and phase change (melting [8,10,22] and boiling [7]).
The phase field method employs a marker function, the phase field φ, to describe the distribution of the various phases. This marker function assumes a constant value in the bulk of each phase (φ = ±1); the interface is represented as a transition layer, across which the phase field changes smoothly between the bulk values. The transport of the phase field variable is described by an advection-diffusion equation; the diffusive term (derived from a Ginzburg-Landau free energy functional) drives the system towards a minimum energy configuration and restores the phase field equilibrium profile across the interface.
The phase field method brings in several important advantages, among which are the implicit tracking of the interface and the automatic handling of topological changes. In addition, the phase field is an Eulerian variable defined in the entire domain, allowing for the adoption of efficient numerical solvers and parallelization techniques. The implicit tracking of the interface decouples the computational cost from the interface extension: with respect to other interface tracking methods, whose computational cost depends on the number of Lagrangian marker points used to track the interface, the phase field method computational cost is independent of the interface morphology and extension. The strengths of the method are mitigated by some side effects that arise when the phobic behaviour (tendency to separate into pure phases) of the two phases is described using a double-well potential in the Ginzburg-Landau free energy functional. The use of a logarithmic or a double-obstacle potential would prevent these effects. However, each of the possible choices has some issues, and the double-well potential is often the preferred one as no singularities, hard to tackle from a numerical point of view, are present.
The adoption of a double-well potential does not bind the diffusion of one phase into the other to the sole interface: each phase can diffuse into the bulk of the other phase (bulk diffusion). Bulk diffusion is the source of shrinkage and coarsening phenomena. Shrinkage occurs whenever the interface profile deviates from its equilibrium (with eventual overshoots/undershoots of the phase field variable) and the volume of the phase enclosed by the interface diffuses into the other phase to restore the equilibrium profile. This way the total mass is still conserved, but there is a mass leakage among the phases. Besides this effect, the diffusive flux introduces also coarsening: since diffusion is driven by an energy minimization criterion, larger domains enclosed by an interface will grow at the expense of smaller ones. In addition to these effects, the velocity field can alter the phase variable interfacial profile; since surface tension forces and thermo-physical properties (density and viscosity) directly depend on the phase field [23,32,33], their accurate representation is affected by alterations of the interfacial phase field profile.
In this paper, we will review and compare different formulations of the phase field method developed to overcome the presented side effects. We consider a classic [19][20][21], a profile-corrected [9,27] and a fluxcorrected [46] formulations (Sect. 2). These are then compared with different benchmarks, including the rise of a bubble in quiescent fluid, the interaction of two droplets in shear flow and the simulation of a large swarm of droplets released in a turbulent channel flow (Sect. 3). Finally, conclusions are drawn (Sect. 4).

Methodology
Turbulent multiphase flows are our focus; therefore, we will present the methods as based on the direct numerical simulation (DNS) of the turbulent field and coupled with the phase field method (PFM) to account for the two phases; in this framework, the Navier-Stokes equation describes the flow field behaviour, whereas the PFM is used to describe interfacial phenomena. The phase field method relies on a marker function, φ, constant in the bulk of the phases (φ = ±1) and undergoing a smooth transition across the interface. In the following, we will briefly present classic [6,19], profile-corrected [9,27] and flux-corrected [46] PFM formulations and their coupling with the Navier-Stokes equations.

Classic phase field method formulation
The time evolution of the phase field φ in the classic formulation is described by a Cahn-Hilliard equation [4,[12][13][14]: where u = (u, v, w) is the velocity field vector, Pe is the Péclet number, and μ is the chemical potential. The Péclet number controls the interface relaxation and is the ratio between the diffusive timescale (proportional to the mobility, here assumed as constant) and the convective timescale. The chemical potential controls the behaviour of the interfacial layer and is defined as the variational derivative of a Ginzburg-Landau free energy functional, F [φ, ∇φ], on a domain : The above expression of F[φ, ∇φ] describes an isothermal immiscible binary mixture and is given by the sum of two different contributions. In particular, f 0 (φ) is a double-well potential, which accounts for the tendency of the system to separate into two pure phases. The second term, Ch 2 2 |∇φ| 2 in Eq. (2), is a non-local term (mixing energy) accounting for the energy stored in the interfacial layer (surface tension). The Cahn number, Ch, is a dimensionless parameter which denotes the thickness of the interfacial layer. Adopting the above-mentioned Ginzburg-Landau free energy functional, the chemical potential results in: At the equilibrium, the chemical potential will be constant throughout all the domain. The equilibrium profile for a flat interface can be obtained solving ∇μ = 0 and leads to: where s is a coordinate normal to the interface. The Cahn number Ch sets the thickness of the interfacial layer where φ undergoes a smooth transition between the two bulk values φ = ±1 following the hyperbolic tangent profile. A mean-field approximation [26] is adopted in the description of the interface: the interfacial thickness (set with the Cahn number) is much larger than that of a real interface (molecular scale). The description of all the scales, from the problem scale down to the molecular scale, would require computational resources currently unavailable and, thus, the adoption of a mean-field approximation. Equations (1) and (4) define the classic PFM formulation [6,19,23].

Profile-corrected phase field method formulation
The profile-corrected PFM formulation [9,27] introduces a penalty flux term, f p , in the right-hand side of the Cahn-Hilliard equation: The penalty flux forces the interface profile towards its equilibrium, reducing the diffusive fluxes induced by the gradients of the chemical potential μ. The penalty flux f p is defined as: with λ being a numerical parameter whose value can be set according to the scaling λ = α/Ch [27], with α being a positive constant. The gradient flow, which includes the penalty flux, is obtained by taking the variational derivative of a free energy functional [27]. Overall, the profile-corrected formulation aims to improve the description of surface forces and thermo-physical properties and to reduce shrinkage and coarsening phenomena, driving the interfacial profile towards its equilibrium.

Flux-corrected phase field method formulation
The flux-corrected formulation [46] aims to further improve the performance of the profile-corrected one. In addition to the penalty flux previously defined, an additional flux, f f , is introduced. This new flux exactly cancels out the component normal to the interface of the flux generated by the chemical potential gradients. This leads to the following Cahn-Hilliard equation: where the new flux introduced is defined as: The flux f f is given by the projection of the chemical potential gradient onto the normal to the interface, which in the frame of the PFM can be computed as n = ∇φ/|∇φ| [5,40].

Coupling with Navier-Stokes equation
In order to describe the flow behaviour, the different PFM formulations have been coupled to continuity and Navier-Stokes (NS) equations. In this work, we consider two phases with different density (ρ 1 = ρ 2 , density ratio γ = ρ 1 /ρ 2 ) but equal dynamic viscosity (η = η 1 = η 2 ). The density is assumed to be a function of the phase field [3,16,23]: The density of the carrier fluid (φ = −1) is taken as reference, and thus it assumes the values: ρ(φ) = 1 in the carrier phase and ρ(φ) = γ in the dispersed phase. Under these assumptions, the continuity and Navier-Stokes equations can be written as follows: where u = (u, v, w) is the velocity field, p is the pressure, g is the gravity unit vector, and the last term of the right-hand side identifies the surface tension forces. Here, surface tension forces are calculated using the Korteweg stress tensor, τ c = |∇φ| 2 I − ∇φ ⊗ ∇φ [25]. The dimensionless parameters are the shear Reynolds number Re τ , ratio of inertial to viscous terms, the Froude number Fr, ratio of inertia to gravity, and the Weber number We, ratio of inertia to surface tension.

Numerical method
The governing equations are solved using a pseudo-spectral method that transforms physical variables into wavenumber space through the use of a Fourier representation along x and y (homogeneous directions) and a Chebyshev representation along z (wall-normal direction). The grid has uniform spacing in the streamwise and spanwise directions, while Chebyshev-Gauss-Lobatto points are used in the wall-normal direction. The Navier-Stokes equations are recasted in a velocity-vorticity formulation: a fourth-order equation for the wall-normal component of the velocity, w, and a second-order equation for the wall-normal component of the vorticity, ω z , are obtained. This approach was adopted by Kim et al. [24] for phases with matched density; since in this work we use phases with different densities (Sect. 3.1) we report in the following the details on the numerical approach we adopted. The left-hand side of Eq. (12) is split into a linear and a nonlinear part: All nonlinear terms are collected in term S, which is dealt with explicitly.
The Navier-Stokes equations thus result in The second-order equation for wall-normal vorticity is obtained by taking the curl of Eq. (15), while the fourth-order equation for the wall-normal velocity is obtained by taking twice the curl of Eq. (15).
In a similar manner, the nonlinear terms of the Cahn-Hilliard equation for the phase field are collected in the term S φ . The nonlinear term S φ for the three formulations is, respectively (classic, profile-corrected, flux-corrected): The Cahn-Hilliard equation thus results in: A splitting technique is adopted to improve the stability of the numerical method [6]; the parameter s is the splitting coefficient.
The equations are time-advanced through an IMEX scheme: an implicit scheme is used for the linear terms (Crank-Nicolson), while an explicit scheme is used for the nonlinear terms (Adams-Bashforth for NS and explicit Euler for the three PFM formulations).
We consider a channel-like geometry bounded by two walls located at z = ±1, at which a no-slip boundary condition is enforced on the velocity field: A zero-flux boundary condition is used for both phase field φ and chemical potential μ at the two walls. The boundary conditions on chemical potential can be reworked in a no-flux condition on the phase field wall-normal third derivative: The no-flux boundary conditions adopted lead to the conservation of the integral of the phase field φ over time [44]: where is the computational domain. The numerical scheme presented above has been implemented in a Fortran 2003 proprietary code and parallelized using a pure MPI approach together with a 2D domain decomposition technique. The Fourier and Chebyshev transforms are performed using the latest FFTW libraries.

Results
To characterize the strengths and weaknesses of the different PFM formulations, their performances have been benchmarked via 2D simulations and then tested considering a 3D simulation of a droplet-laden turbulent flow. The 2D benchmarks consist of a rising bubble in quiescent fluid and of the collision and coalescence of two droplets in a shear flow. Then, starting from the outcome of these first benchmarks, the performances of the formulations have been analysed on a swarm of large and deformable droplets released in a turbulent channel flow.

Bubble rising in a channel
In this section, we consider a bubble that rises in a two-dimensional channel, as sketched in Fig. 1. The bubble, diameter d = 1, is released at the bottom of the channel in an initially quiescent fluid, u = 0. The pressure gradient in the y-direction balances exactly the weight of the fluid (gravity acts in the y-direction), so that the only motion is due to the bubble buoyancy. The density ratio is set to γ = 0.1, while the dynamic viscosity is the same for both phases, η 1 = η 2 = η. The channel size is L y × L z = 2π × 2; the domain is periodic in the y-direction, whereas z is the wall-normal direction. Reynolds, Weber and Froude numbers are set to Re τ = 10, We = 1.0 and Fr = 0.1, respectively. The domain is discretized with N y × N z = 512 × 513 grid points. In order to accurately resolve the thin interfacial transition layer, the Cahn number (dimensionless thickness of the interfacial layer) is set to Ch = 0.02. The Péclet number is set using the scaling Pe = 1/Ch = 50 [28]. The parameter λ, present in the corrected formulations, is set using the scaling: λ = α/Ch, with α being a positive constant. The sensitivity of the formulations to the parameter λ has been tested considering three different values of α, specifically α = 0.02 − 0.05 − 0.1. The bubble, initially at rest, starts to rise along the y-axis. Its rise induces a mean shear that deforms and can eventually break the bubble interface. The time behaviour predicted by the three formulations is reported in Fig. 2 for λ = 5.0 (α = 0.1). We first consider the classic formulation, panels (a, d, g). The bubble, initially circular, deforms, panel (g), and rises along the y-direction, panels (d, a). First, we can notice that, especially for t = 0.5 and t = 0.75, the thickness of the interfacial layer shows some variations along the interface. The layer becomes thinner in the upper part of the bubble and thicker at the trailing edge. It can be also noted that the phase field φ in the bulk exhibits some overshoots (dark red regions). In addition, as the bubble increases its speed and approaches the upper part of the channel, panels (d, a), the two tails become more elongated and separate from the main body of the bubble. This breakage produces a series of smaller bubbles which are then shrank and dissolved in the bulk (light blue region).
This behaviour can be compared against those predicted by the profile-corrected formulation, panels (b, e, h), and by the flux-corrected formulation, panels (c, f, i). The early stages are similar, panels (g, h, i), but after the separation of the tails from the main body of the bubble, t 0.5, some differences can be observed. The smaller droplets are not anymore shrank, panels (b, c); indeed, the additional fluxes introduced in the corrected formulations reduce shrinkage phenomena. The behaviours predicted by the corrected formulations are similar; only few differences can be observed in the later stages of the simulation, panels (b, c). In particular, adopting the flux-corrected formulation, the trailing droplets are slightly bigger and already merged in two larger droplets (bottom part of the tails). This small differences can be addressed to the additional flux present in the flux-corrected formulation, which further reduces shrinkage.
To quantify these qualitative observations, the instantaneous profiles of the phase field φ and the time behaviour of the dispersed phase mass have been computed. In particular, in Fig. 3, the instantaneous profiles of the phase field at t = 0.5 along the y-axis (dash-dotted line of Fig. 1) are shown and compared against the theoretical profile, panels (a-c). We start from the classic formulation, panel (a). At the bubble head, the The shrinkage of the smaller droplets generated by the tail breakage has been quantified computing the dispersed phase mass, m(t), normalized by its initial value, m 0 . In the initial stages of the simulation, for all the formulations, mass leakages are negligible (m(t)/m 0 1). After t 0.5, the tails of the bubble undergo several breakage phenomena and mass leakages can be observed; the mass leakage is larger for the classic formulation (up to 12% at t = 1.0) and largely due to the shrinkage of small bubbles. The fluxes introduced in the corrected formulations are able to suppress shrinkage phenomena: indeed, for the fluxcorrected formulation, at t = 1.0, the mass leakage is about 2.5% (m(t)/m 0 0.975), while it is slightly larger for the profile-corrected formulation, about 4% (m(t)/m 0 0.96).
Lastly, we investigate the effect of the penalty flux considering different values of the constant α and, thus, of λ. We performed two further simulations for each corrected formulation with different values of λ, in particular λ = 2.5 (α = 0.05) and λ = 1.0 (α = 0.02). Overall, three different values of λ were tested: λ = 1.0, 2.5 and 5.0. To show the effect of λ, the instantaneous profiles of the phase field sampled along the centreline of the channel (dash-dotted line in Fig. 1) at t = 0.5 are reported in Fig. 4. To better appreciate the effect of the parameter λ, only the central part of the domain has been reported. The profile obtained from the classic formulation (black solid line) is reported as reference. For the profile-corrected formulation, panel (a), for decreasing λ, the profile shows a slight overshoot in the front of the bubble and φ is less uniform in the bulk. Indeed, for decreasing λ, the penalty flux is reduced and thus similar results as those of the classic formulation are obtained. In the flux-corrected formulation, λ = 2.5 performs best: with decreasing λ a slight overshoot is observed at the bubble front, while with increasing λ small oscillations are observed in the rear part of the bubble.
From the mass leakage point of view, the performances of the corrected formulations deteriorate decreasing λ: for the profile-corrected formulation the mass leakage is larger (5% for λ = 2.5 and 8% for λ = 1.0), whereas for the flux-corrected one the additional flux f f helps in keeping the mass leakages limited (3% for λ = 2.5 and 3.5% for λ = 1.0).  Overall, considering the problem of a rising bubble, the corrected formulations allow for a better conservation of the interfacial profile during the computation. This leads to a reduction in shrinkage phenomena (and thus a better mass conservation of each phase) and to a better representation of surface tension forces and viscosity/density contrasts. Lastly, we would like to remark that this is a challenging benchmark from the mass leakage point of view; indeed, high shear stress values are found along the interface and, as pointed out by Yue et al. [45], they play a crucial role in the shrinkage phenomenon. Moreover, the breakage of the bubble tails generates a large number of small bubbles and further amplifies the mass leakage [45].

Droplet-droplet interaction in shear flow
The accurate description of topological changes is one of the main advantages of the phase field method [31,39,45]. Here, the capabilities of the different formulations in describing topological changes have been compared. Specifically, we consider the coalescence of two droplets in a laminar shear flow. A sketch of the configuration is reported in Fig. 5: two circular droplets, diameter d = 0.7, are placed at y c = π ∓ y/2 and z c = ± z/2, respectively. The shear flow drives the droplets one towards the other; after an initial approaching stage, in which the thin liquid film is drained, the two droplets coalesce.
We assume two phases with matched density (ρ = ρ 1 = ρ 2 ) and viscosity (η = η 1 = η 2 ). The channel has dimensions L y × L z = 2π × 2, and the walls move along the y-axis with opposite directions, v(±1) = ±1. The y-component of the velocity is initialized with a linear velocity profile, v(z) = z; the other velocity components are set to zero, u = w = 0. The Reynolds and the Weber numbers are set to Re τ = 0.5 and We = 1.5; gravity is neglected (no density differences). The domain is discretized with N y × N z = 512 × 513 grid points. The Cahn number is set to Ch = 0.02 and the Péclet number to Pe = 1/Ch = 50 [28]. The parameter λ is set using the scaling λ = α/Ch; even for this case, three values of α, and thus of λ, have been tested. The behaviour predicted by the three PFM formulations is reported in Fig. 6 for λ = 5.0. The droplets interface is identified by the iso-contour φ = 0; the different panels (a-d) refer to different times, specifically t = 2.0, 2.2, 2.4 and 2.6. Different colours identify the formulations: classic in black, profile-corrected in blue and flux-corrected in red.
At t = 2.0, panel (a), the droplets, initially circular, have been deformed by the shear flow and have been shifted one towards the other. All the formulations predict a similar shape; only small differences can be noticed between the classic formulation and the corrected ones. Indeed, some small mass leakages have already occurred and the droplets (classic formulation) are slightly smaller. Later on, at t = 2.2 panel (b), the droplets are closer and the difference among the three formulations becomes clearer: the penalty flux accelerates the coalescence process [46]. As soon as the two thin transition layers overlap, the droplets start to merge. A similar behaviour is observed in the flux-corrected formulation; in this case, since the component normal to the interface of the flux induced by the chemical potential is cancelled, this effect is less pronounced. Then, in panels (c, d), the droplets merge and surface tension forces reshape the droplet. Overall, for λ = 5.0, the corrected formulations slightly anticipate the coalescence. This can be attributed to two main factors: (i) classic formulation exhibits a small mass leakage (up to 2% against almost 0% of the corrected ones), so droplets are slightly smaller and thus the gap between them is wider. (ii) the penalty flux slightly favours the merging and as soon as the two interfacial layers overlap coalescence takes place [46].
The effect of the penalty flux has been investigated considering different values of λ: 5.0 (α = 0.1), 2.5 (α = 0.05) and 1.0 (α = 0.02). In Fig. 7, the time behaviour of the droplet interaction for the profile-corrected, panels (a, b), and for the flux-corrected, panels (c, d), formulations are reported. For the sake of clarity, only the cases λ = 5.0 and λ = 1.0 are shown.
In all the panels, the interface has been identified by the iso-contour φ = 0; black lines refer to the classic formulation (reference case). In panels (a, b), results are marked in blue for λ = 5.0 and in cyan for λ = 1.0, whereas in panels (c, d) results are marked in red for λ = 5.0 and in orange for λ = 1.0. The left column refers to t = 2.2 (before coalescence), whereas the right one refers to t = 2.4 (after coalescence). We start considering the profile-corrected formulation, panels (a, b); we can observe that when λ is decreased (cyan), the difference between the profile-corrected and the classic formulation is reduced. Indeed, reducing the penalty flux magnitude the classic formulation is recovered. As previously noticed, the penalty flux, whose magnitude is proportional to λ, slightly anticipates coalescence [46]. When the flux-corrected formulation is considered, Overall, considering the coalescence between two droplets in shear flow, a fair agreement between the results obtained from the different formulations can be observed. In particular, the penalty flux can influence the coalescence and slightly favours the merging of the two droplets (in particular for λ = 5.0). This effect is less pronounced for the flux-corrected formulation, where the normal flux induced by the chemical potential gradients is cancelled out and, for the lowest value of λ, coalescence is delayed.
From the mass leakage point of view, for all the values of λ considered (corrected formulations), no mass leakages have been observed, whereas for the classic formulation the final mass leakage is around 2%.

Large deformable droplets in wall-bounded turbulence
The results obtained in the previous benchmarks gave us important indications on the performance of the different formulations. Exploiting these indications, three large-scale simulations of a droplet-laden turbulent flow, one for each formulation, were performed.
The formulations have been tested considering a swarm of 256 large and deformable droplets released in a turbulent channel flow. We assume that the dispersed phase (droplets) and the continuous phase (carrier flow) have the same density (ρ = ρ 1 = ρ 2 ) and viscosity (η = η 1 = η 2 ). The droplets, diameter d = 0.4, are released in a channel with dimensions L x × L y × L z = 4π × 2π × 2; the streamwise (x) and the spanwise directions (y) are periodic, whereas z is the wall-normal direction. The initial velocity field is obtained from a direct numerical simulation (DNS) of a fully developed turbulent channel flow at a shear Reynolds number Re τ = 300. The Weber number is set to We = 1.5, the shear Reynolds is kept at the same value (Re τ = 300), and gravity is neglected (droplets and carrier fluid have the same density). The flow is driven by a constant mean pressure gradient in the streamwise direction. The domain is discretized with N x × N y × N z = 512×256×513 grid points; this grid resolution allows for an accurate description of all the flow spatial and temporal scales. The Cahn number is set to Ch = 0.02 and the Péclet number is set to Pe = 1/Ch = 50 [28]. The parameter λ is set from the scaling λ = 0.05/Ch = 2.5 (α = 0.05) for both the corrected formulations. This value was identified from the previous benchmarks, as it allows for an optimal conservation of the dispersed phase mass  and of the interfacial profile (and thus it well represents surface tension forces and density/viscosity contrasts) without significantly affecting coalescence.
A qualitative view of the system is reported in Fig. 8; the droplets interface (cyan) is identified by the iso-contour φ = 0. From Fig. 8, it can be observed how turbulence modifies the shape of the droplets, which are not anymore spherical. The snapshot reported refers to t = 0.5; at this stage, droplets are essentially transported by the carrier flow. Later in time, turbulence can drive droplets towards coalescence or breakage events [33,36].
The effectiveness of the corrected formulations in preserving the interfacial profile during all the computation is of crucial importance. Adopting the classic formulation, the turbulent flow field strongly perturbs the equilibrium profile; the consequent diffusive flux that restores the interfacial profile causes mass leakage between the two phases. By opposite, when the corrected formulations are used, the additional fluxes well preserve the interfacial profile during the entire simulation and the mass leakage is limited. To quantify the effectiveness of the different formulations, the behaviour of the dispersed phase mass, m(t), normalized by its initial value, m 0 , has been computed and reported in Fig. 9. The different lines refer to classic (black), profilecorrected (blue) and flux-corrected (red) formulations. The flux-corrected one exhibits the best performance: the mass leakage is almost negligible (0.5% at the end of the simulation). A similar result is obtained with the profile-corrected formulation: in this case, the mass leakage is slightly higher (up to 2.5% at the of the simulation). It can also be observed that most of the mass leakage for the corrected formulations occurs in the early stage of the simulation, and then the dispersed phase mass keeps constant. As the classic formulation shows a large mass leakage, an additional simulation with halved Cahn number (Ch = 0.01) was performed. A lower Cahn number reduces mass leakage among the two phases, but it must be noted that reducing the interface thickness (the Cahn number) increases the grid requirements. Indeed, in this simulation the grid had to be refined: the number of grid points in each direction was doubled. The dispersed phase mass for this latter case is reported in Fig. 9 with a black dashed line. The lower Cahn number helps in preserving the dispersed phase mass; nonetheless, the corrected formulations exhibit an even lower mass leakage with a higher Cahn number (coarser grid).
Overall, considering a droplet-laden turbulent flow, the corrected formulations improve the description of the interfacial profile, leading to a better representation of surface tension forces, density and viscosity contrasts and a better conservation of the dispersed phase mass.

Conclusions
In this work, the performances of the classic phase field method and of two corrected formulations (profilecorrected [9,27] and flux-corrected [46]) have been tested and compared. These formulations, introducing additional fluxes in the Cahn-Hilliard equation, strongly reduce shrinkage phenomena and thus the mass leakages. The formulations have first been benchmarked considering a rising bubble and the coalescence of two droplets in shear flow. These tests have been used to obtain important indications on the effect that the additional fluxes have in the description of coalescence phenomena, in the conservation of the interfacial profile and in the mass leakage among the phases. Having identified the best choice for the parameter λ, which controls the magnitude of the additional fluxes, the capabilities of the formulations have been compared in the simulation of a swarm of large and deformable droplets released in a turbulent channel flow.
Overall, the corrected formulations allow for a better conservation of the interfacial profile during the computation; this leads to a strong reduction in shrinkage and coarsening phenomena and thus to an improved mass conservation of the dispersed phase. In addition, the representation of surface tension forces and of density/viscosity contrasts is also improved. The corrected formulations proved to be effective even when a turbulent multiphase flow is considered, reducing mass leakages, improving the conservation of the interfacial profile and avoiding the use of a finer grid (thus reducing the required computational cost). From the computational point of view, considering the pseudo-spectral solver adopted, using the corrected formulations increases the computational cost by 6% (profile-corrected) and 11% (flux-corrected), since additional fluxes have to be computed.