A discontinuous Galerkin approximation for a wall–bounded consistent three–component Cahn–Hilliard flow model

We present an efficient high–order discontinuous Galerkin (DG) discretization for the three–phase Cahn–Hilliard system. Among the different approaches, we use the model derived in [Boyer, F., & Lapuerta, C. (2006). Study of a three component Cahn–Hilliard flow model], where the consistency is ensured with an additional term in the chemical free–energy. The model considered in this work includes a wall boundary condition that allows for an arbitrary equilibrium contact angle in three–phase flows. The model is discretized with a high–order discontinuous Galerkin spectral element method that uses the symmetric interior penalty to compute the interface fluxes, and allows for unstructured meshes with curvilinear hexahedral elements. The integration in time uses a first order IMplicit–EXplicit (IMEX) method, such that the associated linear systems are decoupled for the two Cahn–Hilliard equations to be solved. Additionally, the Jacobian matrix is constant, and identical for both equations. This allows an efficient resolution of the two systems by performing only one LU factorization of the size of the two–phase system and two Gauss substitutions. Finally, we test numerically the accuracy of the scheme with a convergence analysis, the captive bubble test and a study of two bubbles in contact with a wall.


Introduction
The study of multiphase flows is of broad interest for both the scientific community and industrial applications (e.g. oil, gas, and nuclear industries). Multiphase flows study the evolution of two or more immiscible fluids, which tend to segregate and be separated by a thin interface. In this work, we study the dynamics of three dissimilar coexisting phases.
The mathematical modeling of multiphase flows distinguish two main approaches: the sharp and the diffuse interface methods. On the one hand, the sharp interface method considers an infinitely thin interface, which is tracked using, for example, a level-set method [1,2]. Then the fluid dynamics equations for each individual phase are solved, being two fluids coupled through an interface condition which enforces mass, momentum, and energy balances [3]. On the other hand, in the diffuse interface (or phase field) methods [4,5] the interface is regularized and it is provided with a finite, yet thin, thickness. As a result, the concentration of the phases, and all the thermodynamic properties, vary smoothly across the interface.
Among the different diffuse interface methods, the Cahn-Hilliard type models are attracting a lot of attention [6,7,8,9,10,11,12,13]. The Cahn-Hilliard equation integrates the effects of phase separation (segregation) and phase homogenization (mixing) in a free-energy function, which is minimized as the flow evolves. The Cahn-Hilliard model for two phase flows is the simpler, and it has been widely studied in the past [14,15,12,13]. Extending the two phase model to three phases is not immediate, since new physical considerations should be taken into account: the model has to be consistent, i.e., if one of the three phases is not present initially, it can not emerge in further times. Besides, it is desirable that the model reduces to the two-phase Cahn-Hilliard model when one of the phases is not present. In this work, we use the model presented in [16], which solves the consistency problem considering a particular choice of the chemical free-energy. For N ⩾ 4 phases, other methods that involve degenerate diffusion coefficients have been developed [10,17,18]. The model is augmented with the boundary condition developed in [19], which permits the prescription of the contact angle between the different phases and the wall by solving a non-homogeneous Neumann boundary condition. The three phase model is numerically approximated in space with a highorder Discontinuous Galerkin Spectral Element Method (DGSEM) [20] that uses the symmetric interior penalty method [21]. The DGSEM is desirable since it allows arbitrary order of accuracy [22,20], the representation of arbitrarily complex geometries through the use of unstructured meshes with curvilinear elements [23], efficient mesh adaptation techniques [24] and the design of provably stable schemes [25,26,27,12,28,13]. The three component Cahn-Hilliard flow model has been discretized by means of the finite element method [29], local discontinuous Galerkin method [30] or spectral element method [10]. The DGSEM has been used by the authors to discretise the two component Cahn-Hilliard flow model [12]. To the authors known this is the first implementation of the three component Cahn-Hilliard flow model in a discontinuous Galerkin Spectral Element method framework.
Finally, we consider a first order IMplicit-EXplicit (IMEX) time integrator. We use an IMEX method since the Cahn-Hilliard equation features a linear fourth order spatial operator (which is solved implicitly), and a nonlinear second order spatial operator (which is solved explicitly). Therefore the solution of the fully-discrete system involves the solution of one linear system for each of the Cahn-Hilliard equations. The two linear systems, however, are decoupled and such that the Jacobian matrices are constant in time and identical for both Cahn-Hilliard equations. This method is efficient as permits a resolution approach where only one LU factorization is performed for the two equations.
The rest of this work is organized as follows: in Sec. 2 we describe the three component Cahn-Hilliard model. The construction of the discrete DG approximation is described in Sec. 3. Lastly, we provide numerical experiments in Sec. 4 that assess the accuracy of the method. Final conclusions and discussions can be found in Sec. 5.

Model description
The phase field approach to multiphase flows introduces one scalar field c j per fluid, which represents the relative concentration of phase j (i.e. the volume occupied by phase j divided by the total volume) in each space-time point. The phases conservation condition, then, reads, In this work, we restrict ourselves to three phase flows, N phases = 3.
We detail in this section the three-phase model derived in [16]. As usual in phase field methods, each of the concentration fields is subjected to a Cahn-Hilliard diffusion equation. The Cahn-Hilliard equation is such that the evo-lution of the concentration minimizes a free-energy. In [16], the free-energy function is, The free-energy contains two terms. The first term contains the chemical freeenergy, F [σ] 0 , which in this model is a polynomial function on the three concentrations c = (c 1 , c 2 , c 3 ). The second term is the interfacial energy, where σ is the interface tension tensor, whose entries are the surface tension coefficients of all the possible interfaces, σ ij , The interface tension coefficients σ ij are positive constants. Finally, ε is the interface width parameter, which controls the thickness of the diffuse interfaces.
The Cahn-Hilliard equation is constructed such that the concentration time derivatives are proportional to the Laplacian of additional scalar fields called chemical potential, µ = (µ 1 , µ 2 , µ 3 ), where the positive constant M 0 is called mobility, and Σ j are constant coefficients called spreading factors. Although in [16] the authors show an extension that allows negative spreading factors, this work only considers positive spreading factors. For the three-phase model, the chemical potentials are defined as, For simplicity, we define f i as to write the chemical potential as µ i = 12 The model is completed with the definition of the chemical free-energy F [σ] 0 . In [16] it is done so that it satisfies two important properties: Property 1 Conservative model. If the initial condition satisfies c 1 (x, 0) + c 2 (x, 0) + c 3 (x, 0) = 1, then, for all further times a conservative model maintains This property is ensured by an appropriate construction of the chemical potentials (5). We sum the three Cahn-Hilliard equations to find that, Therefore, it suffices that the chemical potentials satisfy, which they do by construction. The relation (9) allows us to compute the chemical potential of one of the phases (typically the third phase) as a function of the other two phases, Property 2 Consistent model. If in the initial condition c j (x, 0) = 0, then As shown in [16], this property is guaranteed with the chemical potential construction (5), and an appropriate choice of the chemical free-energy.
We are left with two properties to be satisfied by the chemical free-energy to obtain a consistent model. First, the chemical free-energy of a two-phase system (when one of the phases is not present in the flow, e.g. phase three) is, A consistent chemical free-energy must satisfy that if one of the phases is not present, the resulting chemical free-energy is equivalent to that of a two-phase flow for the other two fluids. A natural approach that satisfies this property is to add the chemical free-energies for the three possible pairs, Second, the reduction to a two-phase chemical free-energy is not enough to guarantee the consistency of the model. If c 3 (x, 0) = 0, then, to ensure c 3,t = 0, the chemical potential µ 3 must be algebraically zero. This implies that the chemical free-energy must satisfy that Eq. (14) is automatically satisfied if an additional term is added to the nonconsistent chemical free-energy (13), The chemical free-energy (15) completes the three-phase model. Note that this additional term also cancels when one of the phases is not present, and the chemical free-energy still reduces to that of a two-phase model (12). As a result, one Cahn-Hilliard equation (4) per phase is solved, with the chemical potential defined in (5) and the chemical free-energy in (15). In practice, since (7) holds, we only solve two Cahn-Hilliard equations (4), and the concentration of the third phase is obtained from the constraint (1), Its associated chemical potential, if needed, is also computed from the other two phases using (10).

Boundary conditions
The Cahn-Hilliard equation is complemented with an appropriate choice of the boundary conditions. Here we adopt a wall boundary condition model with non-zero contact angle for the three-phase Cahn-Hilliard equations.
Since the equation is fourth order in space, two boundary conditions must be specified, for both the chemical potentials and the concentration [4,5]. To guarantee the phases conservation, we enforce a homogeneous Neumann boundary condition for the chemical potential, For two-phase flows, a wall boundary condition with contact angle can be found in [31]. The latter is achieved with a non-homogeneous Neumann boundary condition for the concentration, The contact angle convention is represented in Fig. 1, such that θ ij represents the angle of an interface between fluids i and j with fluid j, and θ ji is that with fluid i. The relation between both angles that was applied in (17) is The stress equilibrium at the wall contact point implies that, For three-phase flows, one can not choose freely the three contact angles with the wall, since a constraint relates them. We write the three wall stress equilibrium equations, Now we subtract the second equation to the first, and then add the third equation to the result, to find that wall contact angles and interface tension coefficients are constrained by the relation Lastly, in [32,11] the non-homogeneous Neumann boundary condition was extended to three phases preserving the consistency property, Eq. (21) allows us to conveniently set the wall equilibrium contact angles.

Discontinuous Galerkin approximation
In this section we present a novel high-order discontinuous Galerkin approximation of the model presented in [16]. Among the different advantages of using a DG method, we highlight the capability to construct numerical approximations with arbitrarily high order of accuracy, the use of curvilinear elements in unstructured meshes to represent complex geometries and the design of provably stable schemes.

Approximation and differential geometry
The physical domain Ω is tessellated with non-overlapping hexahedral elements. Then, we use a polynomial transfinite mapping to geometrically transform the elements to a reference element E = [−1, 1] 3 . The mapping relates coordinates through a function x = X (ξ). The details on how this function is constructed for a general curvilinear hexahedral can be found in [12].
The transfinite mapping is used to transform curvilinear elements to the reference element E. Associated to the transformation, we can compute the covariant and contravariant basis, and the Jacobian, (22) By construction, the contravariant basis satisfies the continuous metric identities, namely, The contravariant basis is also used to transform the differential operators involved in the equations. The gradient and divergence operators are, (see [23]), where ∇ ξ = (∂/∂ξ, ∂/∂η, ∂/∂ζ), and M is the metrics matrix, Additionally, the product of a vector f by the transpose of the metrics matrix M is called a contravariant vector,f .
In the reference element, E, the solution is represented by polynomials of degree N, using tensor product Lagrange interpolating polynomials, l j (ξ), The Lagrange polynomials are written in a set of Gauss or Gauss-Lobatto points, {ξ i } N i=0 . With the definitions (24), we approximate the solution inside the reference element as, where the time-dependent nodal degrees of freedom are The convention used in this work is to adopt lower cases for continuous solutions, and upper cases for the polynomial approximations. The approximation of the metrics is done differently to (27) to obtain a discrete version of the metric identities (23). Following [23], we compute the contravariant basis as, whereê i is the unit vector along the i-th spatial direction.
Based on the Lagrange interpolating polynomials and the choice of the interpolation points, we construct Gauss quadrature rules to approximate integrals inside the reference element. We define the inner product of F and G as their product integral in the reference element, The numerical quadrature weights w ijk = w i w j w k are computed from the one dimensional integrals of the Lagrange polynomials, The numerical quadrature exactly approximates the integral of degree 2N ± 1 polymomials (+1 for Gauss, -1 for Gauss-Lobatto). Therefore, since the Lagrange polymomials satisfy the cardinal property, l i (ξ j ) = δ ij , they are discretely orthonormal in the reference element, where δ ij is the Kronecker delta.
Lastly, given the reference element, the associated surface integral of a contravariant vectorf = ( f ξ , f η , f ζ ) extends to all six faces that define the element, wheren is the outward pointing normal vector to the six faces of the reference element, dS ξ is the surface differential in reference space. The discrete approximation of surface integrals is computed using the quadrature rules in two dimensions, We can represent surface integrals in both physical and reference space. The relation between physical and reference integration variables is, where i = 1, 2, 3 and J i f = |Ja i | is the face Jacobian. Moreover, from the definition of contravariant fluxes we get, Therefore, we can write,

Discontinuous Galerkin approximation of the three-phase model
The Cahn-Hilliard equation is fourth order in space. Thus, to construct a DG approximation, we rewrite the two Cahn-Hilliard equations as a four equation system of first order equations per phase. To do so, we introduce the auxiliary variables g c,i = ∇c i and g µ,i = ∇µ i so that, Recall that f i was defined in (6). Next, we transform the gradient and divergence operators to the reference space using (24), We construct one weak form for each of the four equations. To do so, we multiply by arbitrary N order polynomials, φ i , and integrate in E, Note that in the second and fourth equations we have moved the metrics matrix to the vector test function to obtain contravariant test functions. Finally, we integrate by parts all integrals containing differential operators, and write the resulting surface integrals in physical space using (36), We introduce the polynomial approximation ansatz in (40). We replace continuous functions by the order N polynomials and the exact integrals by the quadrature rules, Since no continuity requirements to the discrete solution have been imposed at the inter-element faces, the boundary integrals are not uniquely defined. Thus, we replace the solution at the inter-element boundaries by an unique solution, represented with the star. The inter-element solution and fluxes are responsible of the coupling between adjacent elements and the enforcement of boundary condition at the physical boundaries. In this work, we use the symmetric interior penalty method (IP), where the gradients are locally computed using (24), and σ is the penalty parameter computed with the estimation by [33], The evolution equation for the coefficients is obtained replacing the test function by the Lagrange polynomials φ = l i (ξ) l j (η)l k (ζ).

Boundary conditions
The boundary conditions enforcement is performed through the numerical fluxes at the physical boundaries. For the chemical potential gradient, we use a homogeneous Neumann boundary condition. Thus, the chemical potential is taken from the interior element, and the normal gradient is set to zero, For the concentration, we apply the non-homogeneous Neumann boundary condition that accounts for arbitrary wall contact angle (21). Thus, we use the interior value for the concentration, and the normal gradients are set to, As in the continuous boundary condition, the wall contact angles θ w ij are subjected to the constitutive constraint (20).

Time discretization
The semi-discrete scheme (41) is complemented with the numerical integration of the left hand side time derivative coefficients. Looking at the continuous equation (4), with chemical potential (5), one finds a linear bi-Laplacian operator for the concentration c i , and a non-linear Laplacian term corresponding to the chemical free-energy derivatives. The bi-Laplacian operator involves fourth order derivatives, which translates in a severe numerical stiffness that restricts the time-step size in explicit solvers. Since the bi-Laplacian operator is linear, but the Laplacian of the chemical potential is non-linear, a commonly adopted technique is to use an IMplicit-EXplicit (IMEX) method [18,12,13].
We revisit the continuous setting (4) to describe the approximation in time. We use the IMEX version of the first order Euler method described in [18]. We define the coefficients c n = c(t n ) as the solution evaluated in t n , such that we evaluate the chemical free-energy in the old time step, t n , and the interfacial energy in the new time step, t n+1 , The stabilizes the non-linear terms, where S 0 is a positive constant, while retaining first order accuracy. Note that the addition of the stabilization term modifies the Jacobian of the implicit solver, but maintains the linearity in c n+1 i . Note that this implementation has two advantages: 1. The implicit system is decoupled: for each of the two phases, an individual linear system only involving c n+1 i is solved. 2. The Jacobian matrix of the implicit solver is identical for the two phases.
Hence, only one Jacobian matrix needs to be computed and stored. This matrix is constant in time due to the linearity.
In this work, the linear system of equations is solved with an LU factorization and Gauss substitution. Since the Jacobian matrix is constant in time, and the LU factorization is done only one time at the beginning of the computations, the resulting implementation is very efficient. However, the algorithm does not restrict to other techniques, e.g. iterative solvers.
Introducing the IMEX scheme (47) in the semi-discrete DG formulation (41), we get the fully-discrete system, In (48), the θ superscript refers to IMEX variables, i.e., those which are not evaluated either in n or in n + 1, but a combination of both (the chemical potential, and its gradient).

Numerical experiments
In this section we perform numerical experiments to evaluate the scheme presented and its numerical implementation. We first study the accuracy of the scheme with a convergence analysis, which solves a manufactured solution. Then, we study the captive bubble test, where a bubble is immersed between the two other phases. Lastly, we test the wall contact angle boundary condition solving two bubbles of two different phases immersed in the third phase.

Convergence analysis
We perform a two-dimensional convergence analysis based on the manufactured solution used in [18] to solve four phase flows, (1 + cos (4πx) sin (4πy) sin (t)) , with the physical parameters given in Table 1. The physical domain is Ω = [−1, 1] 2 , and we enforce periodic boundary conditions at the four physical boundaries. We solve the fully-discrete system (48) in a Cartesian 2 2 grid until a final time t F = 0.1, to then measure the L 2 errors as, The L 2 errors for both fluids concentration are represented in Fig. 2. The  (N > 11), where the error stagnates. We see that the error is higher for c 1 than for c 2 , and that the behavior is not regular, with systematically higher convergence rates for odd polynomial orders. We see that for c 2 , odd polynomial orders yield smaller errors than the theoretical linear behavior, while for c 1 , on the contrary, even polynomial orders yield higher errors than expected. We do not have an answer for this behavior, although it has been seen in other works that particular choices of the manufactured solution might lead to this even-odd phenomena [34,25,12].

Captive bubble simulation
In the second test, we study a bubble (of phase 3) immersed in two layers of the other two phases. As a result of the interface tension, the equilibrium is reached when the angles between the different phases at the triple point (c 1 = c 2 = c 3 = 1/3) satisfy [16], We consider the three different cases for the interface tension coefficients described in [16], which are summarized in Table 2 along with the equilibrium triple point angles computed with (51). For all the three cases, the initial which is represented in Fig. 3(a) condition in all physical boundaries (with θ w ij = 90 • wall contact angles). Recall that phase 1 is the upper layer represented in color white, phase 2 is the lower layer represented in color black, and the bubble is phase 3, colored in gray. The rest of the physical and numerical parameters are given in Table 3. The final equilibrium solution is represented in Fig. 3 for the three cases stud- ied. In the first case ( Fig. 3(b)), given σ 13 < σ 23 , the bubble rises from the initial position, as a result of a higher interfacial tension of the bubble with phase 2 than compared with phase 1. In the other two cases (Figs. 3(c),3(d)), the conditions are symmetric as σ 13 = σ 23 , and we obtain a lenticular shape. We observe that when σ 12 > σ 13 = σ 23 , the lens flattens. These results are in agreement with those in [16]. Moreover, in the right triple point, we have represented with red lines the angles estimated from the theoretical equilibrium configuration (51). We confirm that all three solutions show a good agreement with the theory (51) and with the reference [16].

Wall contact angle simulation
The last test case assesses the exactness of the wall contact angle boundary condition (21). We consider two bubbles of phases 1 and 2 in contact with the inferior and immersed in phase 3. We specify the equilibrium wall contact angles θ w 13 and θ w 23 , and compute the remaining third angle θ w 12 from (20). We maintain the same domain and mesh (with polynomial order N = 8) used for the captive bubble simulation. The new initial condition is represented in Fig. 4(a). The rest of the physical parameters are given in  Table 4. We consider three different cases for the wall contact angles, given in Table 5, and such that each one results in a different configuration of the bubbles in equilibrium. The results are represented in Fig. 4. In red line, we have represented the specified angles, following the convention described in Fig. 1. The first case, Fig. 4(b), represents a configuration where phase 1 is hydrophilic (wets the wall) and phase 2 is hydrophobic (rejects the wall). In the second case, Fig. 4(c), both phases are hydrophobic, whereas in the third case, Fig. 4(c), both phases are hydrophilic. All in all, the contact angles with the wall are in agreement with those imposed by the boundary condition.

Conclusions
In this work, we have developed and implemented an efficient high-order DG scheme for the three-phase Cahn-Hilliard equation. The model developed in [16], complemented with the wall boundary condition developed in [19] (that allows the prescription of the wall contact angle for each of the three phases) has been chosen. The continuous system of equations is approximated with a discontinuous Galerkin spectral element scheme that uses the symmetric interior penalty method to compute the diffusive fluxes. We a first order IMplicit-EXplicit (IMEX) method to integrate in time, such that non-linear terms are solved explicitly, and linear terms are solved implicitly. The solution of the fully-discrete system involves the solution of one linear system for each of the Cahn-Hilliard equations. The two linear systems, however, are decoupled and such that the Jacobian matrices are constant in time and identical for both Cahn-Hilliard equations. This method is efficient as permits a resolution approach where only one Jacobian matrix is needed and only one LU factorization is performed for the two equations. Finally, numerical experiments have been performed to evaluate the scheme presented and its numerical implementation. The introduced method presents spectral convergence and correctly reproduces the contact angles between the different phases (both in the bulk and in the wall).
grant (EUIN2017-88294, Gonzalo Rubio). This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 785549 (FireExtintion: H2020-CS2-CFP06-2017-01). The authors acknowledge the computer resources and technical assistance provided by the Centro de Supercomputación y Visualización de Madrid (CeSViMa).