Journal of Fluid Flow, Heat and Mass Transfer (JFFHMT)

ISSN: 2368-6111

1. Introduction

Phase field methods [1-2] have been successfully used to simulate the flow of two or more immiscible fluids including head-on droplet collision [3-4], droplet impact on solid surface [5-6], and dripping to jetting transition [7-8]. As opposed to sharp interface methods, the phase field approach describes the interface as a transition region of physical quantities with a small but finite thickness, and it is capable of handling complex topological transition naturally without ad-hoc procedures [9]. Of the phase field models for two immiscible, incompressible and density-matched fluids, the basic model so-called Model H [10] has been widely used. It is a system of coupled Cahn-Hilliard and Navier-Stokes equations, and the thermodynamic consistency of this model was proven by Gurtin et al. [11].

There has been a growing interest in the extension of Model H for two fluids where the bulk densities and the bulk viscosities are not matched. Broadly there are two branches of models for such mixture, depending on the definition of mean velocity of two fluids. A quasi-incompressible phase field model for binary mixtures based on the mass-averaged velocity was presented by Antanovskii [12]. Lowengrub et al. [13] extended Antanovskii's model by considering the dependence of the chemical potential on the kinematic fluid pressure. However, with this definition the velocity field is non-solenoidal, and therefore the coupling of the Cahn-Hilliard and Navier-Stokes equations is mathematically much stronger, compared with Model H. For instance the pressure explicitly influences the chemical potential, and the linearization of the two equations becomes much more complicated. For these reasons, there are no discrete schemes available for the full model of Lowengrub, although there is a simplified version of the model used in numerical studies [14-15].

The formulation with the volume-averaged velocity is an alternative to Lowengrub's. This model is based on the assumption of incompressibility, and therefore the velocity field is divergence-free. The incompressibility has not only advantages with respect to numerical simulation, but also maintains the fundamental rule for diffusion, i.e., diffusive flow is controlled by the local compositions, not by the densities, imposed in the Cahn-Hilliard equation, while the total mass of the system is still conserved. Several models with the solenoidal velocity field are proposed with different expressions for the effects of density and viscosity differences (see for example [16-19]), and Aland et al. [20] compared these models for a benchmark example in terms of convergence and accuracy. In addition, various special discretization schemes are presented for these models. The advection term in the Cahn-Hilliard equation was discretized by an upwinding finite volume scheme in [16] and by the finite element method with Muramn and Resetarinera schemes in [17]. Shen et al. [18] compared the formulations with the gauge-Uzawa scheme and the pressure-stabilization scheme for large density ratios.

The objective of this study is to present the numerical scheme with the least-squares spectral element method to solve the incompressible Cahn-Hilliard and Navier-Stokes system for two-phase flow with large density and viscosity ratio. In our previous works we have presented the simulations of the Cahn-Hilliard equation and Model H using the least-squares method [21-22]. The use of the least-squares method has several advantages: (1) it converts the original non-self adjoint operator into a self-adjoin operator, leading to a symmetric positive definite discrete linear system that can be solved by well-established iterative solvers: (2) the inf-sup (or LBB) condition is naturally satisfied, so that equal-order elements can be used [23].

There are few works on two immiscible, incompressible fluids with the least-squares method. Villegas et al. [24] solved two-phase flow equation based on the level-set method with the least-squares method. Tiwari et al. [25] presented an implicit projection method based on the least-squares particle method to deal with two incompressible Navier-Stokes equations for two phases. More recently, Adler et al. [26] used the least-squares method to model Allen-Cahn-type equation.

The fourth-order Cahn-Hilliard equation can be split into two second-order partial differential equations by introducing the chemical potential as an auxiliary variable. In this work, the global differentiability of the set of these two equations and the Navier-Stokes equation is satisfied by approximating the solution with C¹ Hermite polynomial functions in each element. By using higher order global differentiability in local approximations we can avoid introducing excessive auxiliary equations and auxiliary variables and achieve the improved accuracy for the same number of degrees of freedom compared to C⁰ approximation [27]. Surana et al. [28] showed the advantages of the higher order continuity approximation in the least-squares process. For the Cahn-Hilliard equation the use of higher order continuity approximation can result in significant improvement in the accuracy, because the solution rapidly varies over the interfacial region where a thinner interface is preferred to reproduce the correct physics. Discontinuities can result in oscillations in particular when high order polynomials are used, however, it has been demonstrated that the numerical diffusion can be completely eliminated by mesh refinement and p-enrichment in several works [23, 29].

The rest of this paper is organized as follows. In Section 2, we review the governing equations. The numerical approach is described in Section 3. Convergence analysis of our numerical formulation can be found in Section 4 and more general numerical examples are presented in Section 5. We draw conclusions in Section 6.

2. Formulation of Problem

2.1. Notations

The formulations of the phase field method depend on the model of Helmholtz free energy . Cahn et al. [30] proposed a sum of the local Helmholtz energy and surface tension effects on the interface as,

where is the interfacial parameter, and C the concentration, defined as the volume fraction of one phase, i.e., . The local Helmholtz energy has the form of double-well to describe immiscible two phases, and in this study it is approximated as a fourth-order polynomial for isothermal cases [31]:

The total density of the mixture ρ and the viscosity of the mixture μ are defined as

respectively, where and denote the bulk density values and and are the bulk viscosities. Each fluid has its own velocity field, and , and the volume-averaged velocity is defined as

2.2. Governing equations

Following Ding et al. [16], we consider the following coupled Cahn-Hilliard and Navier-Stokes system for two immiscible, incompressible fluid flows with a given density and viscosity ratio:

where M is the mobility, P the pressure, the interfacial tension, and Ω the domain. Here, is the chemical potential for constant temperature, which is the derivative of the free energy with respect to the configuration change, and is the external body force expressed as , with the gravitational acceleration g. In Eq. 6 is the advection term and in Eq. 9 is the capillary force. Eq. 6-9 are along with the following no-penetration and no-slip boundary conditions:

where is the outward vector at the boundary Γ.

In this article, the governing equations are non-dimensionalized. The dimensionless density and viscosity for a liquid-vapor mixture are defined as

where the density ratio and viscosity ratio are defined as the ratios of vapor properties over the liquid one as, and , respectively. Other dimensionless variables are:

where L₀, t₀, U₀, P₀, and M₀ are the reference length, time, velocity, pressure, and mobility. Dropping out the primes, Eq. 6-9 in dimensionless form can be written as

As for the dimensionless groups, Pe is the Peclet number, Re is the Reynolds number, Ca is the capillary number, and Bo is the Bond number defined as

Note that the dimensionless numbers regarding to hydrodynamics are based on the liquid values to see the effect of liquid properties.

The study of the influence of mobility is out of the scope in this study, but two mobilities are used depending on the problems: and . The first mobility is not degenerate in the pure phases, but the second one is degenerate, i.e., zero mobility in the pure phases. Due to this difference, with the polynomial mobility the interface motion is determined by diffusion along the interface itself, i.e., surface diffusion, while the interface motion with the unity mobility is also influenced by a diffusion process in the two bulk phases [32]. The polynomial mobility has been introduced by Cahn et al. [33] from thermodynamic grounds.

Equilibrium of a droplet surrounded by quiescent vapor requires the curvature of the interface to be uniform. Otherwise, the pressure and the capillary forces at the interface are unbalanced, giving rise to a nonzero velocity field. The equilibrium interface profile can be found by minimizing with respect to , i.e., zero chemical potential. With the double-well free energy, expressed as Eq. 2, the equilibrium interface profile can be determined analytically as,

with the z-coordinate chosen along the gradient of C.

3. Numerical Methods

3.1. Least-squares method

The basic idea of the least-squares method is the minimization of the residual functional in a least-squares manner. A set of partial differential equations can be represented as

with is a linear partial differential operator, is a boundary operator, U is a solution vector, and and are corresponding source terms. We assume that the system is well-posed and the operators and are continuous mappings from the solution in Hilbert space onto the data in Hilbert spaces . The least-squares functional corresponding to Eq. 22 can be defined as the square of the residual as,

The solution is calculated from the following minimization statement based on variational analysis:

Equivalently, it is possible to write the necessary condition as:

with

where is a symmetric, positive definite bilinear form, and is a continuous linear form.

3.2. Spectral element discretization

The computational domain Ω is divided into Ne non-overlapped sub-domains , called spectral elements. The local solution in each element is approximated as Eq. 28 by the linear combination of a set of continuous basis functions

with the local coordinate of in the parent element, the unit cube , with d = dimΩ, and the coefficients in the approximation. The global approximation in Ω is constructed by connecting the local solutions, i.e.,

The discretization is based on a space-time coupled formulation with the time-stepping procedure suggested by Pontaza et al. [34]. The solution is approximated on consecutively aligned space-time strips domains, and a strip is composed of only one element in time, with the time step size . The initial condition of the system is applied to the first strip, and for the subsequent strips, the final solution from the previous strip is prescribed to the next strip as the initial condition.

3.3. Iterative method for coupling and nonlinear terms

Eq. 16-19 are the coupled Cahn-Hilliard and Navier-Stokes equations with nonlinear terms. To handle the complexities of nonlinearity and coupling, two iterative schemes are used. The Newton linearization method is used to cope with nonlinear terms, and a relaxation method is used in decoupling of the equations to assure convergence and acceleration of the iterations.

Through the Newton linearization method the nonlinear terms in Eq. 17 and Eq. 19 can be rewritten as the linear forms for the two-dimensional spatial domain, respectively:

The terms with subscript l represent the values from previous iterative steps for linearization, and they are regarded as constants in the linear form.

With the linearized forms Eq. 30-31, the operators and solution vectors in Eq. 22 for the Cahn-Hilliard and Navier-Stokes equations are given in Eq. 32 and Eq. 33, respectively. The terms from previous decoupling step are denoted as subscript n.

Both nonlinear and decoupling convergences are declared when the relative norm of the residual, i.e., , is less than 10^-6, with the residual defined as

3.4. C¹¹¹ approximation

In this work, for a two-dimensional space and time domain , equal order C¹ p-version hierarchical approximation functions, C¹¹¹ approximation, are used to interpolate the solution . The solution space is spanned by polynomial functions of order and which are continuous and differentiable for two space coordinates and the time coordinate, and their first derivatives are continuous:

For constructing the three-dimensional basis functions, the p-version hierarchical interpolation basis functions and construction approach by Solin and Segeth [35] are used. The same basis functions and construction approach have been used in [36] for the solution of the Cahn-Hilliard equation in one-dimensional space and time domain.

The one-dimensional basis functions consist of a set of four vertex basis functions defined as the cubic Hermite polynomials. The corresponding functions in the reference domain [-1,1] are defined as

Among these four functions, and contain the information of the value of the approximation, and and contain the information of the derivative at the boundary of the element. These basis functions can be extended in a hierarchical manner for :

Since these shape functions vanish on the boundary of the element, they are called bubble functions.

The three-dimensional basis functions for the space and time domains are written as the tensor product of one-dimensional basis functions, i.e. . To simplify the numerical manipulation, we can define the three-dimensional basis functions as

with where . Thus, we have basis functions .

3.5. Fully discrete model

The assembling of the element matrix can be written as

The integral expression is numerically approximated by the Gaussian quadrature based on the GLL-roots, and over-integration with larger number of quadrature points Q than the expansion order p is used to improve the convergence rate of the solution [37], .

We implement the conjugated gradient method with the Jacobi preconditioner to solve the system of algebraic equations, Eq. 22, in element level.

This element-by-element technique is suitable to solve a linear system possessing a symmetric positive definite matrix and well-matched with iterative methods, such as for linearization and the time-stepping procedure [38]. A Matlab code developed at our group has been used as a main setup. For parallelization of the algorithm Matlab MPI is used to allocate sub-domains to processors and communicate between processors. The governing equations in each sub-domain are solved separately, and the solution vectors are assembled in a main processor.

4. Code Verification

4.1. Convergence analysis based on the manufactured solution

In this section the presented least-squares spectral element formulation is verified with the manufactured solutions in terms of error indicator. The following manufactured solution, expressed as products of trigonometrical functions, is used:

The manufactured solution is not the analytic solution, so the residuals have to be added to the right-hand side of the discretized equations. Initial and boundary conditions are given by the manufactured solution; Dirichlet boundary conditions for the velocity fields and Neumann boundary conditions for the phase field variables are applied. In the convergence analysis, L²-norm of difference between the approximated solution and the manufactured solution , i.e., is chosen as error indicator and it can be written as

The used parameters are: , and . The governing equations are solved over a domain , with until 3.0 s (6 periods). The time step size is determined through the Courant-Friedrich-Lewy (CFL) condition

where the stable time step size is generated when CFL is around unity. A CFL spanning from 0.05 to 1.6 is used in this study with U₀ of 0.5.

The errors in time for and are presented in Figure 1. Errors oscillate with the period of solution, and after the first period the oscillation amplitudes become stable. It illustrates the fact that the errors do not diverge, but just fluctuate depending on the contour of solutions.

The p-refinement study ranging from to 10 and with and h-refinement study of and with are conducted, and the errors at 3.0 s are presented in Figure 2. The L²-errors in the p-refinement study show exponential decay with respect to the expansion order p, as expected for the least-squares method. The h-refinement study shows the expected linear convergence rate with slope 5, as theoretically predicted for a fixed order .

4.2. Performance of parallel computation

In this work the parallel computation using the element-by-element technique is implemented for a problem with a refined grid. The performance of the parallelization is measured with the manufactured solution example. Two parameters, the speedup ratio and the parallel efficiency defined as Eq. 50 and Eq. 51, are used. The speedup ratio is the ratio between the elapsed time of one processor and the elapsed time of np processors , and the parallel efficiency is the speedup ratio over the number of processors np.

Figure 3 presents the behavior of the speedup and parallel efficiency with 1, 2, 4, 8, 16, 32, 48, and 64 processors for and with and elements and 4 for the expansion order. For all simulations 32 GB of RAM were used. For all cases the parallelization is less efficient as more processors are used because the time required for communication becomes more dominant. Based on these results, the number of 16 processors was adopted for further simulations in this work.

5. Falling Droplet

We model a falling droplet under gravity force in a vertical channel to investigate the effects of density ratio and viscosity ratio of the two fluids. Three parameters are used to show the effects quantitatively; center of mass , circularity c, and falling velocity , which are defined as

with the diameter of projected area . Circularity is defined as a ratio the perimeter of projected-area-equivalent circle to the perimeter of droplet. It is assumed that the thin interface between the phases is spatially located at .

The spatial domain is chosen as a vertical channel, , and the no-penetration and no-slip boundary conditions are applied to all boundaries. Initially a droplet with diameter D = 0.5 is placed at under zero flow field, and the initial interface condition follows the analytical solution, Eq. 21. To assure the equilibrium state numerically as well, the simulation runs without any external force, until L^¥-norm of difference of concentration in time is lower than , i.e., , and after then the gravity force is imposed and the characteristics of the falling droplet are measured.

We set , and for all simulations presented in this work. The simulation runs until the droplet reaches with the time step size , and the mesh size is taken to be with a grid discretization of and , corresponding to CFL of 0.48. The solution is approximated by basis functions with an expansion order 4. An example of a falling droplet with and is presented in Figure 4. We observe that the droplet falls and turns into an ellipse in time, while the interface width and the area of the droplet are preserved.

5.1. Effects of density ratio

Figure 5 shows the variations of (a) center of mass and (b) falling velocity of the droplets with four different density ratios: , and 0.001, while the bulk vapor density is fixed at 1 and the viscosity ratio is fixed at 0.1. The plateau on the falling velocities of and 0.05 cases represent that the droplets have reached their terminal velocities, which are around 2.6 and 3.3, respectively. On the other hand, the falling velocities with and 0.001 show rather linearly increasing velocities whose slopes are 9.0 and 9.6, respectively, which are close to the value of gravity force, and it implies that the drag force for these cases is negligible compared with buoyancy. Due to the linearity of the falling velocity for cases and 0.001, the center of mass in time are well-fitted with the quadratic formulas, and , respectively, while those of and 0.05 show linear regions as the droplets reach their terminal velocities.

The shapes of droplets with different around are presented in Figure 6(a), and their circularities in time are given in Figure 6(b). We observe there are insignificant differences in the droplet shapes with respect to . The buoyancy on the droplet is approximately proportional to when vapor density is much smaller than liquid density. On the other hand, the drag force onto the droplet is proportional to because the same circularity for different represents the same projected area. Therefore it can be concluded that the force balance between buoyancy and drag forces is mainly dependent on , and whether the droplet reaches the terminal velocity in a short time or not is also determined by .

5.2. Effects of viscosity ratio

Figure 7 presents the variations of the (a) center of mass and (b) falling velocity of the droplets with four different viscosity ratios: , 0.5, 0.1, and 0.01, while the bulk vapor viscosity is fixed at 1 and the density ratio is fixed at 0.1. The droplets with and 0.5 show their peak velocities at and 0.5 s, respectively, while in the other two cases the falling velocities approach their peaks. The measurements of center of mass are in accordance with the results of falling velocity by showing the quadratic region at the beginning and a linear region at the end. The center of mass graph for and 0.5 cases after t = 0.4 s are enlarged and presented in Figure 7(a), and it illustrates that the slope for has a slightly concave shape, corresponding to its decreasing velocity after the peak.

Figure 8(a) shows the snapshot of four droplets with different when they pass around , and the deformation of droplets by the location are presented in Figure 8(b), and their circularities in time are recorded in Figure 8(c). The circular shape of droplets in and 0.01 cases are relatively conserved, while in the other two cases the droplets are stretched horizontally and even a dimple at the top is developed. With increasing , the bulk liquid viscosity is lower, and it leads to decrease of the relative strength of the surface tension to the inertial forces. Hence, the interface can endure less pressure difference and the droplet has a more pronounced deformation.

With respect to the correlation of velocity and droplet shape, we can conclude that as long as the falling droplet maintains its elliptic shape, as the cases of and 0.01, its velocity still approaches the peak, but above this point, as the droplet is spread out or forms a dimple, as the cases of and 0.5, the velocity starts to drastically decrease. This difference on the evolution of velocities can be explained by the force balance, as similar to the previous analysis. In this example the increased drag force by primarily results from larger project area as well as higher drag coefficient, as reported in [39-42], by considering the falling velocity and are almost the same for different . The spread shape deaccelerates the droplet, as if it plays a role as the parachute of a skydiver.

5.3. Flow around falling droplet

Figure 9 shows the velocity field around the droplet with and in the right-half domain of the axis of symmetry, . As the droplet moves a doughnut-shaped ring forms behind the droplet (from ), and this ring is separated from its tail in the form of vortex loops. Van Dyke [43] also showed from the experiment that the vortex starts to form for Re larger than 130 and scales of the detached eddy increase considerably with further increase in Re. Since then, the droplet takes a shape of ellipsoid and it leads to growth of circulating region at the rear of the droplet. With the enhanced vortex size the velocity field stretches horizontally the fluid elements on the droplet surface. Through this positive feedback, the droplet is further deformed. The same result about the positive feedback was obtained in [43].

6. Conclusion

We presented the least-squares spectral element scheme for two immiscible, incompressible fluids with large density and viscosity ratio. Two-dimensional coupled Cahn-Hilliard and Navier-Stokes equation with the volume-averaged velocity was solved, and C¹ cubic Hermite polynomials were used in approximating the solution. A time-stepping procedure, Newton linearization and element-by-element technique were used. The p- and h-refinement study with the manufactured solution was conducted to verify our solver, and the performance of parallelization was measured. A falling droplet was numerically handled with our solver - the effects of density and viscosity ratio on the droplet dynamics were studied, and the flow around the droplet was analyzed. According to the classification of [44], all the simulation cases in this study are sorted as 'sheet thinning' state, which is expected to have breakup with long time travel. Therefore, for future study, studies on longer simulations or higher viscosity ratio to observe sheet-like shape of droplet and ultimately break up into small pieces are recommended.

Acknowledgements

This work is supported by the Research Council of Norway (FRINATEK project 21442), Norway.

References