An efficient phase-field method for turbulent multiphase flows

With the aim of efficiently simulating three-dimensional multiphase turbulent flows with a phase-field method, we propose a new discretization scheme for the biharmonic term (the 4th-order derivative term) of the Cahn-Hilliard equation. This novel scheme can significantly reduce the computational cost while retaining the same accuracy as the original procedure. Our phase-field method is built on top of a direct numerical simulation solver, named AFiD (www.afid.eu) and open-sourced by our research group. It relies on a pencil distributed parallel strategy and a FFT-based Poisson solver. To deal with large density ratios between the two phases, a pressure split method [1] has been applied to the Poisson solver. To further reduce computational costs, we implement a multiple-resolution algorithm which decouples the discretizations for the Navier-Stokes equations and the scalar equation: while a stretched wall-resolving grid is used for the Navier-Stokes equations, for the Cahn-Hilliard equation we use a fine uniform mesh. The present method shows excellent computational performance for large-scale computation: on meshes up to 8 billion nodes and 3072 CPU cores, a multiphase flow needs only slightly less than 1.5 times the CPU time of the single-phase flow solver on the same grid. The present method is validated by comparing the results to previous studies for the cases of drop deformation in shear flow, including the convergence test with mesh refinement, and breakup of a rising buoyant bubble with density ratio up to 1000. Finally, we simulate the breakup of a big drop and the coalescence of O(10^3) drops in turbulent Rayleigh-B\'enard convection at a Rayleigh number of $10^8$, observing good agreement with theoretical results.

on the same grid. The present method is validated by comparing the results

Introduction
Turbulent multiphase flows are ubiquitous in nature and technology. Examples are raindrops [2,3], ocean waves [4], fuel sprays [5], and the transmission of virus-laden droplets during respiratory events [6,7,8], just to name a few. In order to gain deeper insights into their complex and rich behavior, efficient, high-fidelity computations are crucial. For turbulent multiphase flows, direct numerical simulations (DNSs) present far greater challenges than for single-phase flows [9]. The reasons are the much finer length-scales and faster time-scales induced by the existence of the second phase, especially when the deformable interfaces between the fluids break up or coalesce.
In the phase-field method, two immiscible phases are represented by their volume fractions C and 1 − C, respectively. The spatial distribution of C is determined by the Cahn-Hilliard equation [35,20,36]: The quantity in square brackets is the chemical potential defined by the variation of free energy with respect to C. It includes an excess free energy term (the first term), and a bulk energy term (the second term) with ψ = 1/4 C 2 (C−1) 2 being the simplest non-singular form that has two equal energy minima, namely at C = 0 and 1 [35,20,36]. Physically, ψ represents the bulk energy density due to the inhomogeneous distribution of volume fraction in the interfacial region. We will give more technical details in Section 2.1. In Eq. (1), the Laplacian of the first term on the right hand side is biharmonic, i.e. it contains fourth-order derivatives. State-of-the-art solvers for the standard single phase flow Navier-Stokes equation is highly efficient and well-studied for turbulent flows. This is because the typical algorithm to solve the computationally demanding Poisson equation-a necessary step for enforcing incompressibility-is based on fast Fourier transforms (FFTs) [37,38], as described in Ref. [39]. In Refs. [40] and [1], the FFTs is extended to multiphase flows by employing a split method, meaning the variable-coefficient pressure-gradient term is spilt into an implicit constant term and an explicit variable term. As a result, the Poisson equation can be solved up to 40 times faster than without the split method [1].
However, with the application of FFTs in multiphase flows, the computational cost of the biharmonic term becomes the new bottleneck for the phase-field method. The reason for this is that the common solution technique for the biharmonic term in the phase-field method involves an implicit solution that requires 25 grid points for a second-order spatial discretization, see Fig. 1(a) (details in Section 3.1). Therefore, in this study, we will in particular focus on an optimal discretization of the biharmonic term. We propose a novel discretization scheme for the biharmonic term in the phasefield method to couple with the approximate-factorization method, which is an efficient way to implicitly solve the hyperbolic systems [41] and easily parallelize it. We will implement the phase-field method [36] with this novel scheme into our open-source DNS package AFiD (www.afid.eu) [42,39], which is a second-order finite difference solver that has been well-validated in many studies of turbulent flows [43,44,45]. AFiD is highly-parallelized with a pencil distributed strategy [39,46], and includes an FFT-based Poisson solver [42]. In addition, we will apply a split method [40,1] to the pressure solver to deal with large density differences between the two phases.
To validate the present approach, we simulated cases of drop deformation in a shear flow and of a rising buoyant bubble. Our results are compared to previous studies and are further assessed using a grid convergence study. Finally, we simulated the case of a breakup of one big drop as well as the coalescence of O(10 3 ) drops in turbulent Rayleigh-Bénard convection, and show the good performance of the present approach for large-scale computation.
The paper is organized as follows. The governing equations are introduced in Section 2. Then we address the numerical methodology in Section 3. In Section 4, we simulate several test cases to validate our approach and show its ability to deal with turbulent multiphase flows in large-scale computation. We conclude our study in Section 5.

Cahn-Hilliard (CH) equation
Turbulent flows with two incompressible immiscible fluids are investigated here. We use the phase-field method [47,36] to capture the interface between two fluids. Here, the sharp interface is modeled by a diffused one with finite thickness, and represented by contours of the volume fraction C of fluid 1, and thus the volume fraction of fluid 2 is 1 − C. The evolution of the volume fraction C is governed by the Cahn-Hilliard equation, where u is the flow velocity. We choose the Péclet number (the ratio of advection and diffusion) and the Cahn number (a dimensionless measure of the thickness of diffuse interface) the same as in Ref. [22], i.e. Pe = 0.9Cn and Cn = 0.75h/L with h and L the uniform mesh size and the characteristic length, respectively. To enforce mass conservation, the correction method proposed by [48] is used. This correction method resembles that of Ref. [49] and exhibits good performance (see Section 4.3.1).

Navier-Stokes (NS) equations
The fluid motion is governed by the momentum and continuity equations, which have been made dimensionless using the properties of fluid 1. Here, u is the velocity and P the pressure. ρ and µ are the density and the dynamic viscosity, respectively, which are both functions of C defined as, where λ ρ = ρ 2 /ρ 1 and λ µ = µ 2 /µ 1 are the ratio of the densities and viscosities of the two phases (denoted by the subscript), respectively. The surface tension force F st is computed as in [36], In Eq. (3), the gravity force is G = −ρ/Fr j with j being the vertical direction. The dimensionless numbers controlling the problems are thus the Reynolds number Re = ρ 1 UL/µ 1 , the Weber number We = ρ 1 U 2 L/σ, and the Froude number Fr = U 2 /(gL), where σ the surface tension coefficient, g the gravity acceleration, and U is the characteristic velocity.

Numerical method
We use staggered meshes and solve the CH equation on the uniform mesh with size h for all three directions and the NS equations on the stretched mesh: the procedure for the coupling of the two meshes (uniform and stretched) is based on that reported in [50] and it is described in Section 3.4. A lowstorage third-order Runge-Kutta method [51] is used to temporally advanced all the equations. The biharmonic term in Eq. (2), viscosity term in Eq. (3), and diffusion term in Eq. (27) are implicitly solved by the Crank-Nicolson scheme, while the other terms are solved explicitly. In spatial discretization, central second-order accurate finite-difference schemes are used for all terms (details can be found in [42,36]), except for two: one is the advection term of volume fraction C in CH equation (2), which is solved by fifth-order WENO scheme [36], and the other is the biharmonic term which is solved by a novel scheme proposed in Section 3.1.

Discretization of biharmonic term in CH equation
To accurately advance the CH equation (2) with a large time step, we should implicitly solve the biharmonic term Cn 2 ∇ 4 C at the right-hand side of Eq. (2). At the same time, its discretization scheme should retain the same order of error as the term ∇ 2 (C 3 − 1.5C 2 + 0.5C), which is also at the righthand side of Eq. (2) and discretized by central second-order finite-difference schemes of O(h 2 /L 2 ).
Typically, the biharmonic term is discretized according to Fig. 1(b) (we restrict the expression to a 2D case for the ease of representation), (8) When we implicitly solve this expression, the presence of mixed partial derivatives poses challenges for computational cost and code parallelisation.
To circumvent the use of mixed partial derivatives when solving Eq. (8), we propose a new discretization scheme, which is shown in Eqs. (13), (14) and (15). Thus, we can split this discretization into two one-dimensional parts A x C with C i+m,j and A y C with C i,j+n , which means that only the points on the axes remain (Fig. 1b). Then, we can use the approximate-factorization method (described at the end of this section) to efficiently solve Cn 2 ∇ 4 C implicitly. Our main idea is replacing C i±1,j±1 in Eq. (8) with C i+m,j and C i,j+n (Fig. 1b), where m and n = −2, −1, 0, 1, 2. The replacement is justified based on the Taylor series expansions, where we define C ′ e = (∂C/∂e) i,j with e = x, y, s and τ , so do C ′′ e and C ′′′ e . The directions x and y are the perpendicular axis directions in Cartesian coordinates, and the directions s and τ are obtained by rotating x and y by 45 • . Since the Laplacian operator is rotational invariant, we have so we have the relations, where the first and third-order derivatives are eliminated since the points are symmetrical about (i, j). Thus, C i±1,j±1 can be replaced by C i+m,j and C i,j+n as shown in Fig. 1 Substituting Eq. (12) into Eq. (8), we get the new discretization scheme, where the error O(h 2 /L 2 ) is of the same order as the term ∇ 2 (C 3 − 1.5C 2 + 0.5C) at the right-hand side of Eq. (2). Comparing Eq. (13) and Eq. (9), we get the following pentadiagonal matrix, for 2D, where the values in the first and last rows are determined by boundary conditions. Now, with the convenient form of Eq. (14), the approximatefactorization method can be employed to solve the biharmonic term implicitly. The same idea can be directly extended to three-dimensions, and the points used in the mixed partial derivatives are replaced as shown in Fig. 1(a). Thus, we get the operators, for 3D.
With A x , A y and A z , we can use the approximate-factorization method [41,42] to efficiently solve the following equation with the known q l from the previous time step and unknown q l+1 for the next time step, where E represents the terms calculated explicitly, β is the constant coefficient, A x , A y and A z are discretization operators, and (q l+1 + q l )/2 originates from the Crank-Nicolson scheme. Eq. (16) can be rewritten as, Then we factorize the operators on the left, (18) After factorization, the computation only requires inversions of separate tridiagonal matrices rather than the inversion of a large sparse matrix, which leads to a significant reduction in computation cost and memory [41,42]. Then, Eq. (17) can be solved by the following steps, where the superscript * represents the intermediate parameter. In Eqs. (19), (20) and (21), the inversion of matrix will be extremely cheap when we carefully choose A x , A y and A z , respectively, provided they only involve the points in one dimension.

FFT-based solver with a split method for Poisson equation with large density contrast
The NS equation (3) is solved here by a projection method, where u * is an intermediate velocity field calculated from Eq. (3) without the pressure term. Considering ∇ · u l+1 = 0, we have, To solve this Poisson equation with large density variations, we use the split method proposed by [1] to apply fast Poisson solver to Eq. (23). In the split method [1], the Poisson equation (23) with the variable coefficient 1/ρ l+1 is split into an implicit constant density part and an explicit variable part, where we define ρ 2 ≤ ρ 1 . Substitute Eq. (24) into Eq. (23), Then, a standard fast Poisson solver can be used here. After getting P l+1 , the velocity field is updated as,

Pencil distributed parallel strategy
The parallel method in the present approach is a pencil distributed parallel strategy (details in [39]). Here, the computational domain is split in two dimensions and this strategy allows us to use more CPU cores for large-scale computation, such as 70 billion points with 64K cores as reported in [39]. The other advantage is that this strategy is well coupled with the approximatefactorization method to implicitly solve the equations. The high performance of this parallel method has been extensively validated in [39] and [46]. Moreover, it has already been used in many studies of turbulent flows in large-scale simulations [52,43,44,53].

Multi-resolution meshes for C and u
One feature of our method is that the volume fraction field C can be integrated on a refined uniform mesh, even if the momentum field u is integrated on a non-uniform mesh. For the C field, a uniform mesh is a recommended choice. The reasons for this are as follows: The computation of surface tension force is key to simulate multiphase flows. To ensure the truncation error of surface tension in space is of the same order in all directions, uniform mesh spacing in each direction is necessary near the interface. Furthermore, considering the drops in turbulent flows is likely to break up into smaller sized drops and distribute throughout the domain, the use of a uniform mesh can easily handle the spatially dispersed drops. Therefore, the uniform mesh is a good choice for C field.
On the other hand, in wall-bounded turbulence, the resolution requirements of the u field are more restrictive at the walls, where very thin kinematic boundary layers need to be resolved. The same strict requirements apply for Rayleigh-Bénard convection, where a large number of near-wall nodes are required to resolve thin thermal boundary layers [54]. Therefore, a stretched non-uniform mesh is a good choice for resolving u or the temperature field. This multi-resolution treatment of the mesh allows for large computational savings [50,55] since the operations are by far cheaper when integrating the momentum field on coarser meshes, as compared to the single scalar Cahn-Hillard equation without any elliptic equation.
The multi-resolution method that decouples u and C works as follows. u is projected from a base mesh, which is non-uniform, to a refined uniform mesh on which C resides. The projection employs a tri-cubic Hermite spline interpolation, with a stencil of four points in each direction, for a total of sixty-four points in three dimensions. Here, the Hermitian interpolation is a preferred since the accuracy has been proven to be sufficient for turbulent flows, and is considerably cheaper than other methods such as B-splines [50]. This stencil is generated only once at the start of the simulation and is reused throughout. To preserve the solenoidal properties of the momentum field, instead of directly projecting u, the normal velocity gradients on the base mesh are first computed and then the projection is applied on the normal velocity gradients. Finally, with a refined 2D velocity field interpolated at a reference location (in each direction), the refined velocities are integrated for the entire domain using the interpolated gradients. For the back-coupling of the C field, the refined uniform mesh is directly projected to the stretched mesh since there is no solenoidal requirement for C. This down-sampling projection step is used to obtain µ, ρ and F st . The present method is an improvement over the previous method used in [50], since here, the stretched mesh can contain an arbitrary number of nodes employing different stretching parameters.

Results and discussion
In Section 4.1, we test the convergence of the results with mesh refinement and the performance of the new discretization scheme for the biharmonic term. Section 4.2 shows the ability of the present approach to deal with large density and viscosity contrasts. In Section 4.3, a possible application of multiphase turbulence is simulated -Rayleigh-Bénard convection with drops, where the performance of the multi-resolution meshes is also tested.

Drop deformation in shear flow
In order to test the mesh refinement convergence of our approach and the performance of the new discretization scheme for the biharmonic term, we consider the deformation of a drop in a shear flow with matched density and viscosity. A drop of radius R is initially placed at the center of a domain of 8R × 8R × 8R, as shown in Fig. 2. In the domain, there are two no-slip plates moving at a speed of U in opposite direction, and periodic boundary conditions are used in the other directions. Due to the shear stress exerted by the surrounding fluid, the drop elongates until the surface tension counteracts the resulting load. We define the deformation ratio Γ = (L−B)/(L+B) as in [56,57,58] to quantify the degree of drop deformation, where B and L are the lengths of the minor and major axes of the deformed drop at equilibrium, respectively, see Fig. 2. The governing dimensionless parameters are the capillary number Ca = µγR/σ, the Reynolds number Re = ργR 2 /µ, and the Weber number We = ρ(γR) 2 R/σ = Ca Re, whereγ = 2U/H is the shear rate and H the thickness of the fluid layer. Gravity is not considered here. With Ca ≪ 1 and Re ≪ 1, Γ is expected to linearly depend on Ca accounting to Γ ≈ (35/32)Ca [59].    [59] and the previous numerical results [56] gives good agreement. With increasing Ca, the deformation ratio Γ becomes larger than the theoretical prediction [59] since the assumption of Ca ≪ 1 for this prediction is no longer satisfied. As reported in the previous studies [56,57,58], the drop breaks up at Ca = 0.39 and Re = 1. We also perform this case in a domain of 12R × 8R × 8R as shown in Fig. 4. The drop breaks up into three smaller ones as expected. Fig. 5 shows the results of the convergence study with different mesh size h = 0.0031, 0.0042, 0.0050, 0.0063 and 0.0100 at Ca = 0.1. The numerical error E h is calculated by comparing Γ to the value obtained with the finest mesh (h = 0.0031). The convergence rate is of 1.4, which is between 1 and 2, as expected since the phase-field method [47,36] for the interface used here is first-order accurate while the NS solver is second order [42].
We have also tested the performance of the explicit discretization scheme in Eq. (8) and the new implicit scheme in Eq. (15) for the biharmonic term Cn 2 ∇ 4 C described in Section 3.1. Since the explicit scheme requires a small time step, here we consider the quantity Γ max , which is reached around t = 0.4, instead of Γ at equilibrium, which is attained only around t = 30. Here we show a convergence study with mesh refinement at δt = 5 × 10 −5 (the largest value to maintain numerical stability) with the explicit scheme and δt = 2 × 10 −3 with the new scheme in Fig. 5, where the results agree well. It shows the new implicit scheme in Eq. (15) is highly efficient and accurately discretizes the biharmonic term. Thanks to this, we can perform the largescale simulations of turbulent multiphase flows.

Rising bubble with buoyancy
In this subsection, we test the performance of the present approach by simulating a three-dimensional bubble rising in liquid water with a large density and viscosity contrast up to 1000 and 100 times, respectively, which has the same configuration as previous axisymmetric studies [36,60]. Initially, we place a bubble (fluid 2) of radius R in the domain of 8R×8R×8R with the distance from the bottom plate to bubble center of 1.6R, as shown in Fig. 6. Noslip and non-penetration boundary conditions are enforced at all boundaries. The dimensionless parameters controlling this problem are the Reynolds number Re = ρUR/µ = 100, the Bond number Bo = ρUR 2 /σ = 200, and the density and viscosity ratios λ ρ = 0.001 and λ µ = 0.01, respectively. Note that Fr = 1 and We = Bo due to the characteristic velocity U = √ gR. The mesh used here is 400 × 400 × 400, where the mesh size is the same as in the axisymmetric simulations [36,60].
Thanks to buoyancy, the bubble rises. For Bo ≫ 1, with the surface tension not large enough to counteract buoyancy, the bottom of the bubble will rise faster than the top, as shown in Fig. 7. Therefore, eventually, the bubble breaks up from the tip and evolves into a toroid. Although here a three-dimensional case is performed to test the performance of our code, the flow is indeed axisymmetric, so that we can compare our results with the previous studies of axisymmetric simulations [36,60]. In our numerical simulations, the breakup occurs at t = 1.61 and y = 4.1R, which agrees well with the simulations from previous studies using different numerical approaches, which are t = 1.60 and y = 4.05R with level set method [60], and t = 1.61 and y = 4.09R with diffuse-interface method [36]. Besides, Fig. 7 presents the comparisons of the bubble shape at different time instants t = 0.8, 1.6 and 2.4. It shows that also the shape of the bubble's interface in the present study is in good agreement with the previous ones [36,60].

Multiphase turbulent Rayleigh-Bénard convection
Here we consider a possible application of multiphase turbulent flows by using the present approach: Turbulent Rayleigh-Bénard convection with drops, as shown in Fig. 8. Rayleigh-Bénard convection is the motion of a fluid in a cell heated from below and cooled from above [61,62,63].
For Rayleigh-Bénard convection, the temperature advection equation reads where c p = k d /(κρ) is the specific heat capacity. The thermal conductivity k d is defined as where λ k d = k d2 /k d1 is the ratio of the thermal conductivity. We choose the distance between the hot and cold plates as the characteristic length, and the free fall velocity U = √ α 1 gL∆ as the characteristic velocity. The relevant dimensionless groups of the configuration are the Rayleigh number Ra = α 1 ρ 1 gL 3 ∆/(µ 1 κ 1 ), the Prandtl number Pr = µ 1 /(ρ 1 κ 1 ), where α is the thermal expansion coefficient, ∆ the temperature difference and κ the thermal diffusivity, in addition to the dimensionless numbers controlling the droplets.
For this case the gravity force G in Eq. 3 depends on both C and the dimensionless temperature θ, whose effects on density are considered within the Boussinesq approximation, where λ α = α 2 /α 1 is the ratio of the thermal expansion coefficients α.

Breakup of one big drop in turbulent Rayleigh-Bénard convection
Initially, a drop of radius 0.23H (represented by C = 1) with matched density and viscosity with the ambient fluid is placed at the center of the domain H × H × H, with a linear temperature profile and zero velocity. The boundary conditions at the top and bottom plates are set as C = 0, no-slip condition and fixed temperature θ = 0 (top) and 1 (bottom). Periodic boundary conditions are used in the horizontal directions. The chosen dimensionless parameters are Ra = 10 8 , Pr = 1 and We = 8000. Note that We here is large because it is defined by the system height instead of the droplet size. For local Weber number which is defined using the droplet size, we find that the value is O(1) after the droplet breakup, which is consistent with the Kolmogorov-Hinze theory [64,65]. The chosen Rayleigh number is large enough for the flow to enter the turbulent regime. The mesh is 500 × 500 × 500, which is consistent with the grid resolution checks in [66]. Fig. 9 shows snapshots of the drops in Rayleigh-Bénard convection. The drops first deform due to buoyancy (see Fig. 9a), and then breaks up because of the small surface tension (see Fig. 9b). As time evolves, hundreds of drops of various sizes are advected in the turbulent field ( Fig. 9c and 9d). The drop size is characterized by an effective diameter D, which is defined as   Fig. 10. We observe that the probability distribution function (PDF) of the large drops follows the scaling (D/H) −10/3 while that of the small drops obeys the scaling (D/H) −4/3 , which both originate from the previous theory studies for the respective regimes [67,68]: First, in turbulent flows, the distribution of the drop size has been studied extensively. The well-known −10/3 scaling law for the large drops in turbulence was proposed in ref. [67] and validated by many experimental and numerical studies [69,70,13,17]. Second, the derivation of the −4/3 scaling law for the relatively small drops originates from a recent study [68]. It is based on the energy balance in a regime dominated by surface tension. Fig. 10 shows that the present numerical simulations and the theoretical analyses [67,68] give consistent results.
The mass conservation is also tested in this section. Fig. 11 shows the normalized mass loss l loss = (m t − m 0 )/m 0 , where m t is the mass of fluid 1 (drops) at time t and m 0 is the initial mass of the drop. We see that the maximal mass loss is of the order of 10 −5 and the value of l loss is not increasing in time. This demonstrates the good mass conservation in the present approach, which is consistent with the other studies with phase-field methods [36,71,49].
We also simulated the case on multi-resolution meshes with otherwise unchanged parameters, uniform mesh of 500 3 for the CH equation and stretched  The empty symbols are the present data, and the filled symbols the data of turbulent single-phase flows [46]. We also test the computational efficiency of the method on the supercomputer MareNostrum at the Barcelona Computing Center (2 sockets Intel Xeon Platinum 8160 CPU with 24 cores each @ 2.10GHz, for a total of 48 cores per node). Two sets of gridpoints are used, i.e. 1000 3 and 2000 3 , and the option of multi-resolution is not used here to fit the setting of the previous study. The wall clock time per step and the speedup comparing with a single core as functions of CPU cores are presented in Fig. 12. Compared to the AFiD code for single phase flows [46], the computational cost of the present approach for the multiphase flows is only less than 1.5 times more. Moreover, the parallel efficiency is quite good until the CPU cores used are more than 3072. These data show that the computational performance of the present approach for turbulent multiphase flows is nearly as good as the solver for turbulent single-phase flows.

Coalescence of O(10 3 ) drops in Rayleigh-Bénard convection
The topological change of the interface includes the breakup and coalescence of drops. In Section 4.3.1, we clearly observed the breakup of drops. In this section, we will show the coalescence of O(10 3 ) drops in turbulent Rayleigh-Bénard convection. The initial setup is presented in Fig. 13 we placed 1014 drops with a uniform diameter of 0.08H in a domain of 2H ×2H ×H. The simulation was performed on the mesh of 1000×1000×500 and 2048 CPU cores. The Weber number was set to We = 1000, which is smaller than that in Section 4.3.1. The other dimensionless parameters and boundary conditions are the same as in Section 4.3.1.
As seen from the snapshots at t = 10, 40 and 150 in Fig. 14 (a), most of drops coalesce into larger ones. Since the Weber number here is smaller than that in Section 4.3.1, surface tension here is stronger and can resist inertia, leading to larger drop sizes.
We also simulated a case with a different initialization, where only one big drop with a diameter of 0.8H is placed at the center of the domain. Although different initial conditions are used, similar statistic equilibrium states were obtained after sufficiently long times (see Fig. 14).

Conclusion
In this study we have shown how to efficiently implement the phasefield method into the single-phase DNS solver AFiD. A new discretization scheme for the biharmonic term Cn 2 ∇ 4 C of the Cahn-Hilliard equation has been proposed. Together with the approximate-factorization method, the FFT-based Poisson solver, and a pencil distributed parallel strategy, massive DNSs (up to 8 billion gridpoints and 3072 CPU cores are used) for turbulent multiphase flows can be performed.
The suggested new approach has then been validated by comparisons with several numerical experiments. In the case of drop deformation in shear flow, the results agree well with theoretical and previous numerical results, and the convergence study with mesh refinement shows an accuracy between first and second order, as expected. Then, also for the case of a rising bubble with buoyancy, good agreement is achieved when comparing our results with previous simulations, even with large density or viscosity contrast of up to 1000 or 100 times, respectively. Furthermore, in the case of breakup and coalescence of drops in turbulent Rayleigh-Bénard convection, we observe good performance of our approach to deal with turbulent multiphase flows, including good mass conservation and high efficiency of computation, thus establishing our scheme to perform reliable simulations for turbulent multiphase flows in large-scale computations.
The new scheme and code therefore offer great opportunities to better understand the physics of turbulent two-phase flow with coalescence and breakup of droplets and bubbles.