{'fraction_non_alphanumeric': 0.047339283696757425, 'fraction_numerical': 0.020280194328324484, 'mean_word_length': 4.159720198192947, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A divergence-free HDG scheme for the Cahn-Hilliard phase-field model for two-phase incompressible flow

We construct a divergence-free HDG scheme for the Cahn-Hilliard-Navier-Stokes phase field model. The scheme is robust in the convection-dominated regime, produce a globally divergence-free velocity approximation, and can be efficiently implemented via static condensation. Two numerical benchmark problems, namely the rising bubble problem, and the Rayleigh-Taylor instability problem are used to show the good performance of the proposed scheme.


Introduction
In the past few decades there has been tremendous progress in the development and analysis of numerical methods for single-phase incompressible flow problems. Recently, the importance of producing an exactly divergence-free velocity approximation for single-phase incompressible Stokes and Navier-Stokes equations within the mixed finite element framework has been stressed [9]. There are active research on the extension of numerical methods developed for single-phase incompressible flow problems to two-phase incompressible flow problems [6]. The fundamental issue relevant for the simulation of two-phase flows that is non-existent in one-phase incompressible flow problems is the numerical treatment of the unknown interface.
The function W (φ) is a double well potential, here we use The function M (φ) is a mobility function, here we use where γ is a constant mobility coefficient. Finally, defines a length scale over which the interface is smeared out, and the parameter σ is a scaled surface tension, which is related to the physical surface tension σ by σ := 3 2 √ 2 σ. It is noteworthy to straighten that the model (1), unfortunately, does not admit an energy law. Following [7], an energy estimate of the model could be shown by introducing the new variable r = √ ρ, see [15]. We also mention that thermodynamically consistent phase-field models that admit an energy law was recently derived in [1]. It was numerically shown in [2] that very similar results were obtained for standard finite element discretizations for the model (1) and the thermodynamically consistent model in [1] in the context of a rising bubble benchmark problem [8]. For this reason, we focus on the construction of HDG scheme for the slightly simplier model (1), although it does not admit an energy law.
The rest of the paper is organized as follows. In Section 2, we first introduce the divergence-free HDG-based spatial discretization for the model problem (1), then apply a standard Crank-Nilson based IMEX time discretization for the resulting ODE system, which leads to a (conditionally stable) linear decoupled fully discrete scheme. Then in Section 3, we first numerically study the convergence property of the proposed scheme, which indicates second order convergence in time, and optimal (k + 1)-th order of convergence in space when polynomials of degree k is used for the spatial discretization. Next, we apply our scheme to two benchmark tests, namely, the rising bubble problem, and the RayleighTaylor instability problem. We conclude in Section 4 with some future work.

The divergence-free HDG scheme
In this section, we introduce the divergence-free HDG scheme for the model problem (1) in two dimensions. Although our scheme can be defined on hybrid triangular/quadrilateral meshes, in this paper we only present the scheme on structured rectangular meshes. To this end, let T h := {T } be a conforming rectangular triangulation of the rectangular domain Ω, and let E h = {F } be the collection of edges of T h . We set h to be the maximum mesh size of T h . Given a rectangular element T , we denote P m,n (T ) as the space of polynomials of degree at most m in the first argument and at most n in the second argument on the element T . We further denote Q m (T ) := P m,m (T ) to simplify notation. Given an edge F , we denote P m (F ) as the space of polynomials of degree at most m on the edge F . On each element T , we denote the tangential component of a vector v on an edge F by tang(v) := v − (v · n)n, where n is the unit normal vector on F .
The following finite element spaces will be used: where k is the polynomial degree, and the curl operator in the velocity space (5d) is the rotated gradient: ∇× = (∂ y , −∂ x ). Note that functions in X k h are defined only on the mesh skeleton and are continuous on the mesh nodes, while functions in M k h are defined only on the mesh skeleton and have normal component zero. Note also that functions in the velocity space (5d) is globally divergence-free, and {∇ × ξ j } N j=1 forms a set of basis for V k h , whenever {ξ j } N j=1 forms a set of basis for the scalar continuous finite element space Φ k+1 h . In particular, integration-by-parts yields the following identity: which will help us to exclude the presure approximation in the numerical scheme.
The proposed spatial discretization is given as follows: , and ν h := ν(φ h ) to shorten the notation. Here (·, ·) denotes the L 2 -inner product on the mesh T h , and the operators in (7) are given as follows: , and α > 0 is a sufficiently large stabilization constant to ensure positivity of the viscous operators D h ((c, ·, ·), (·, ·)) and B h ((c, ·, ·), (·, ·)), which is taken to be α = 4 in all our numerical simulations.
We remark that in the above scheme (7), the embeded discontinuous Galerkin (EDG) approach [4,12] is used to discretize the second-order terms in equations (7a) and (7b). The convection term in (7a) is discretized using an EDG-based upwinding. The convection term in the momentum equation (7c) is discretized using a classical DG-based upwinding, the viscous term therein is discretized using a divergence-free HDG approach with projected jumps [10,11] to further save computational cost, the surface tension force term is obtained by integration by parts, and the pressure term cancels out due to the equality (6). In particular, we emphasis that: (i) the numerical dissipation (from upwinding) in the treatment of the two convection terms in (7a) and (7c) is beneficial to stabilize the scheme in the convection-dominated regime without using any extra residual-based stabilization; (ii) due to the choice of the (globally) divergence-free velocity space V k h , the divergence-free condition (1d) is always satisfied strongly for the scheme (7), which is excepted to be more robust than methods that only satisfy the divergence-free condition in an approximate/discrete sense, c.f. [9]. (iii) the introduction of the hybrid unknonws (in the context of EDG for phase-field variables, and HDG with projected jumps for the velocity variable) is a natural way to further improve the efficiency of the resulting linear system solvers (via static condensation). For the temporal discretization, we simply use the second-order Crank-Nilson-Adams-Bashforth IMEX approach [3] as follows. For any positive integer and set t j+1 := t j + δt j , where CF L is the CFL constant. (ii) Extrapolate velocity and phase variables at time t j+1/2 := t j + 1 2 δt j from data at time t j and t j−1 : (iii) Solve the phase field variables at the next time level using the extrapolated velocity: and is a linearization around φ j h . Note that here the unknows φ j+1 h and φ j+1 h stay at time t j+1 , while the unknows µ j+1/2 h and µ j+1/2 h stay at time t j+1/2 . The scheme (8a) is linear, and can be efficiently implemented via static condensation such that the globally coupled degrees of freedom (DOFs) are those on the mesh skeleton only (2 DOFs per vertex, and 2(k − 1) DOFs per dege), whose computational cost is similar to a standard continuous Galerkin finite element method. (iv) Solve the velocity variables at time level t j+1 : In the actual implementaion, we solve for the stream function ξ h ∈ Φ k+1 h and recover velocity via the formula u h = ∇ × ξ h . The scheme (8c) is can also be efficiently implemented via static condensation such that the globally coupled degrees of freedom are those on the mesh skeleton only (1 DOF per vertex, k DOFs per edge for the space Φ k+1 h , and k DOFs per edge for the space M k−1 h . In the practical implementation, we further modify the slave variables ρ and ν as follows: which ensure that

Numerical results
In this section, we present several numerical results for the divergence-free HDG scheme (8) proposed in Section 2. The NGSolve software [13] is used for the simulation.
3.1. Accuracy test. We use the method of manufactory solutions to test the spatial and temporal accuracy of the scheme. In particular, we take the computational domain to be a periodic unit square, add a source term f φ in the equation (1a), and take the source terms f φ in (1a) and f in (1c) such that the exact solution is given as follows: u(t, x) = 0.2 sin(πt) sin(2πx) cos(2πy) −0.2 sin(πt) cos(2πx) sin(2πy) .
Furthermore, we take ρ 1 = 100, ρ 2 = 10, ν 1 = 10, ν = 1, σ = 10, = 0.04, and γ = 10 −3 . The final time is taken to be T = 0.5, and we use a uniform time step size δt = h (k+1)/2 . The history of convergence of the L 2 -norm errors in φ h , µ h , and u h at the final time T = 0.5 on a sequence of uniform square meshes is recorded in Table 1 for k = 1 and k = 2. We clearly observe optimal convergence order of k + 1 in all the variables for both cases, except for the error in µ h when k = 1, where we loose half order of convergence. In particular, this indicates the expected second-order accuracy in time of the temporal discretization.

Rising bubble.
We consider the rising bubble benchmark problem proposed in [8]. The test setup is extensively described in [8]. The domain Ω = [0, 1] × [0, 2] is filled with fluid 1 (φ ≈ 1) except for a circular bubble, which consists of fluid 2 (φ ≈ −1). The initial bubble has a radius of 0.25 with its center at (0.5, 0.5). See Figure 1 for the sketch of the domain and boundary conditions.  Table 2 lists the fluid and physical parameters for the two test cases used in [8]. Because the density of the bubble is smaller than the density of the surrounding fluid (ρ 2 < ρ 1 ), the bubble rises. The evolution of the bubble is tracked for three time units during which the defined benchmark quantities are measured. The measured benchmark quantities are the following: • Center of mass where y is the vertical coordinate. • Circularity: c = perimeter of area-equivalent circle perimeter of bubble = 2 φ<0 π dx φ=0 1 ds .
• Rise velocity: For both test cases, the viscous effects and surface tension forces dominates the convection. We take polynomial degree k = 2 and use a uniform time step size throughout. Uniform square mesh with mesh size h = 2 −5 , h = 2 −6 , and h = 2 −7 are used in the simulation. We take time step size to be δt = 0.005 (600 total steps) for the cases h = 2 −5 and h = 2 −6 , and δt = 0.0025 (1200 total steps) for the case h = 2 −7 . Following [2], we further take the diffuse interface width = 0.64h, and mobility coefficient γ = 10 −3 .

3.2.1.
Results for test case 1. The bubble, being initially circular, first stretches horizontally and develops a dimple at the bottom before it reaches a stable ellipsoidal shape, see Figure 2.  Figure 3 shows the bubble shapes (contour line of φ h = 0) at the final time (t = 3), along with reference data from [8]. It is clear that as mesh size h (and ) decreases, the bubble shape converges to the reference value.   Table 3. Minimum circularity and maximum rise velocity, with corresponding incidence times and final position of the center of mass for test case 1.
Furthermore, we plot the circularity, center of mass, and rise velocity over time in Figure 4. All the quantities seem to converge as mesh size h decreases.

3.2.2.
Results for test case 2. In test case 2, the decrease in surface tension causes the bubble to develop a more non-convex shape and thin filaments. This is a much harder problem, and also in [8], agreement between the used numerical approaches could not be achieved. It even remains unclear if break off of the thin filaments should occur for this setting. We refer to Figure 5 for a depiction of the temporal evolution of the bubble shape at mesh size h = 2 −7 . The initial circular shape is shown to gradually develop two filaments on its sides as it experiences an upward pushing force. But the filaments does not break. This observation is similar to that for the diffusive interface method used in [2]. Figure 6 show the bubble shapes at the final time (t = 3). One can see that, the filaments become thinner for smaller h (and smaller ). Therefore, it is possible that break off happens in the case → 0.
Similar to test case 1, we list in Table 4 the quantitative comparison with the benchmark values. Here, we restricted the comparison with the time interval [0, 2] because also the reference solutions  do not agree well for later times. One can again see that each quantity approaches the reference value as h decreases.  Table 4. Minimum circularity and maximum rise velocity, with corresponding incidence times and final position of the center of mass for test case 2.
Furthermore, we plot the circularity, center of mass, and rise velocity over time in Figure 7. However, even for the finest grid h = 2 −7 , = 0.005, differences remain visible. We conclude that the resolution is still not fine enough to get sufficiently close to the reference solution. 3.3. RayleighTaylor instability. The RayleighTaylor instability is a two-phase instability which occurs whenever two fluids of different density are accelerated against each other. We consider a similar setting as in [7]. This problem consists of two layers of fluid initially at rest in the gravity field in the domain Ω = [0, 1/2] × [−2, 2]. The initial position of the perturbed interface is η(x) = −0.1 cos(2πx). The heavy fluid is above and the density ratio is 3 (ρ 1 = 3, ρ 2 = 1). The (initial) transition between the two fluids is regularized by hyperbolic tangent: We take gravity g = 2 so that the non-dimentional time scale is the same as the time scale of Tryggvason [16]. The viscosity in both fluids is taken to be ν = √ 2/Re, where Re is the Reyholds number. Both Re = 1000, and Re = 5000 are used in the numerical simulations. For the phasefield model parameters, we use = 1.28h, γ = 10 −3 , and take a small surface tension constant σ = 0.01 . We notice that no surface tension effect is taken into account in the model [7]. The same flow boundary condition as the rising bubble problem is used here, namely, the upper and lower boundaries are set to no-slip, and the left and right boundaries are set to slip conditions. Again, we take polynomial degree k = 2 and consider uniform rectangular meshes with mesh size h = 2 −6 and h = 2 −7 . We take variable time step size to be where v max is an estimated maximum velocity magnitude at the current time step. We run the simulation till time t = 2.5. On the coarse mesh with h = 2 −6 , a total of 1808 time steps is used when Re = 1000, and 1966 steps when Re = 5000. On the fine mesh with h = 2 −7 , a total of 3688 time steps is used when Re = 1000, and 4240 steps when Re = 5000.
The time evolution of the contour of the phase-field variable φ h for Re = 1000 and Re = 5000 are plotted in Figure 8 and Figure 9, respectively, at times 1, 1.5, 1.75, 2, 2.25, 2.5. From the results for Re = 1000 in Figure 8, we observe that the solutions on the two meshes are consistent in that they show similar structures and differ only in fine details at large time. Second, by comparing the solutions for the results for Re = 5000 in Figure 9, we observe that the solution on the two meshes are in very good agreement in the early stage of the time evolution t ≤ 1.75. Some noticeable differences occur at later times and consist in the development of structures within the main vortex that are more complex on the fine mesh than on the coarse one. All these results are qualitatively similar to the results in [7], where a projection FEM was used to solve the variable density flow without surface tension.
Finally in Table ??, we compare the minimal and maximum of y position (bubble and spike location, respectively) of the interface φ h = 0 with results in [7] (data in [7] is extrapolated up to two digits accuracy from Figure 4 therein). Excellent agreement of our results with the reference data is observed. In particular, we observe that, on the same mesh, the bubble/spike locations are very similar for Re = 1000 and Re = 5000.

Conclusion and future work
We have presented a novel divergence-free HDG scheme for a Cahn-Hilliard phase-field model for two-phase incompressible flow. The (linear and decoupled) fully discrete is observed to be secondorder accurate in time, and optimal order accurate in space ((k + 1)-th order for the L 2 -errors when polynomials of degree k is used). Benchmark results are presented for the classical bubblerising problem and the Rayleigh-Taylor instability problem, which are consistent with results in the literature.
This work consists of our initial investigation of divergence-free HDG schemes for phase-field model of two-phase incompressible flow. The scheme is robust in the convection-dominated regime, produce a globally divergence-free velocity approximation, and can be efficiently implemented via   Table 5. Bubble and spike locations for the Rayleigh-Taylor instability problem at different time.
we shall consider models admit energy law such as those in [1,15]. The combination of our (spatial) divergence-free HDG scheme with the recently introduced scalar auxilary variable (SAV) technique [14] for the phase-field temporal discretization seems to be a very promising approach, which consists of our ongoing work.