A variational interface-preserving and conservative phase-field method for the surface tension effect in two-phase flows
aa r X i v : . [ phy s i c s . c o m p - ph ] J u l A variational interface-preserving and conservative phase-fieldmethod for the surface tension effect in two-phase flows
Xiaoyu Mao a , Vaibhav Joshi a , Rajeev Jaiman a, ∗ a Department of Mechanical Engineering, University of British Columbia, Vancouver, Canada
Abstract
We present a finite element based variational interface-preserving and conservative phase-fieldformulation for the modeling of incompressible two-phase flows with surface tension dynam-ics. The preservation of the hyperbolic tangent interface profile of the convective Allen-Cahnphase-field formulation relies on a novel time-dependent mobility model. The mobility coeffi-cient is adjusted adaptively as a function of gradients of the velocity and the order parameter inthe diffuse interface region in such a way that the free energy minimization properly opposesthe convective distortion. The ratio of the convective distortion to the free energy minimizationis termed as the convective distortion parameter, which characterizes the deviation of the diffuseinterface profile from the hyperbolic tangent shape due to the convection effect. In the phase-field formulation, the mass conservation is achieved by enforcing a Lagrange multiplier withboth temporal and spatial dependence on the phase-field function. We integrate the interface-preserving and conservative phase-field formulation with the incompressible Navier-Stokesequations and the continuum surface tension force model for the simulation of incompressibletwo-phase flows. A positivity preserving scheme designed for the boundedness and stabilityof the solution is employed for the variational discretization using unstructured finite elements.We examine the convergence and accuracy of the Allen-Cahn phase-field solver through ageneric one-dimensional bistable diffusion-reaction system in a stretching flow. We quantifyand systematically assess the relative interface thickness error and the relative surface tensionforce error with respect to the convective distortion parameter. Two- and three-dimensional ris-ing bubble cases are further simulated to examine the effectiveness of the proposed model onthe volume-preserving mean curvature flow and the interface-preserving capability. Finally, wedemonstrate the applicability of the proposed model for a complex case of two bubbles risingand merging with a free surface, which includes complex topological changes and the surfacetension dynamics using unstructured finite elements.
Keywords:
Phase-field method, Interface-preserving, Conservative, Surface tension,Two-phase flow, Finite elements
1. Introduction
Two-phase flow of immiscible fluids is ubiquitous in many natural phenomena and en-gineering applications. Examples include bubbly cavitating flows around marine propellers ∗ Corresponding author
Email address: [email protected] (Rajeev Jaiman)
Preprint submitted to Journal of L A TEX Templates August 7, 2020 in the numerical simulation of immiscible two-phase flows. When the surface tension forceplays a significant role, the handling of mutual dependency between the representation and theevolution of interfaces becomes considerably challenging. The representation of the interfacegeometry has a direct impact on the surface tension dynamics, which is one of the importantdriving forces during the interface evolution. The quality of two-phase flow solutions is very sensitive to the approximation of surface tension force in the capillary dominated regime. Themodeling of immiscible two-phase flows poses other well-known difficulties regarding highdensity and viscosity ratios, the mass conservation, the discontinuity of properties across theinterface, and the topological changes during simulating realistic flows. In the present work,we particularly focus on the accurate representation of interfaces for the surface tension force calculation in incompressible two-phase flow problems, while retaining the mass conservationduring complex topological changes.The interface between immiscible two-phase fluids can be represented by a sharp or diffuseinterface approach. In the sharp interface approach, the fluid-fluid interface is treated as a sharpboundary separating the domains of the two phases. The boundary is explicitly parameterized by specifying its location and geometry. The parameterization can be accomplished by track-ing the interface with a conforming mesh, which is the approach of the arbitrary Lagrangian-Eulerian method [4] and the front-tracking method [5]. However, the mesh operation duringcomplex topological changes, such as the merging and the breaking-up of the interface, posessignificant difficulties. An alternative approach for the parameterization is to reconstruct the sharp interface according to a function representing the volume fraction [6], which is employedin the volume of fluid (VOF) method [7]. Although the mesh operation is avoided, the volumefraction function suffers from difficulties in the calculation of the normal and curvature fromthe reconstructed interface. This issue is resolved in the diffuse-interface approach by elimi-nating the requirement of interface parameterization and by introducing a smooth phase-field function. In the diffuse-interface approach, the change of physical properties at the interfaceis considered as a gradual variation within a transitional region with a finite thickness. Thesmooth variation can be modeled by a continuous scalar-valued function serving as a phaseindicator for the two-phase mixture. The function is chosen as the signed distance functionto the interface d in the level set method [8], or a hyperbolic tangent function tanh( d/ √ ε ) in the phase-field method, where ε is an interface thickness parameter controlling the thick-ness of the diffuse interface region. With the decrease in the interface thickness parameter,the phase-field function converges to the Heaviside function description in the sharp-interfaceapproach. The gradients of the functions acquire non-zero values only in the diffuse interfaceregion, which can be utilized to reformulate the surface tension force as a volumetric source term. Two popular models to impose the surface tension force are: (i) continuum surface force(CSF) model, which distributes the sharply defined force with a Dirac Delta function to the dif-fuse interface region [9], (ii) free energy based surface tension force model, which imposes theenergy balance in the context of the phase-field method [10, 11]. In this paper, we consider thephase-field method based on the transient Allen-Cahn equation concerning the computational efficiency and simplicity. The CSF model is employed for the surface tension dynamics tocircumvent the chemical potential calculation in the free energy based surface tension model,which is not required in the Allen-Cahn phase-field equation.The hyperbolic tangent shape and the thickness of the interface in the phase-field method211] or in the conservative level set method [12, 13] will not necessarily remain the same during the interface evolution. The diffuse interface region may undergo stretching from the convec-tion, which leads to the development of nonuniformity. We refer to this phenomenon as theconvective distortion. While the thinning of the interface causes the numerical difficulty of re-solving a high gradient with less computational elements, the thickening of the interface resultsin the loss of accuracy. Furthermore, the surface tension force models are subjected to inaccu- racies due to the convective distortion, since their prerequisite interface profiles are no longermaintained. In the level set method, this issue is resolved by solving a reinitialization equationuntil the interface profiles are recovered. This additional procedure can be computationallyexpensive than solving the convection equation [14]. The reinitialization process may also leadto poor mass conservation, which entails further corrections in the level set method [12, 15]. While the phase-field methods resemble the conservative level-set methods, there are some fun-damental differences. For example, the property of maintaining the interface profile close to thehyperbolic tangent function, which minimizes the free energy, is embedded in the phase-fieldmethod. This liberates the phase-field method from the interface reinitialization required in thelevel set method, thereby reducing the computational cost [16] and providing robustness for any kind of geometric manipulation to compute the interfacial curvature. Likewise the level-set method, the phase-field methods have the advantage to handle any topological changes inthe interface due to their Eulerian description. Furthermore, the mass conservation can be en-forced in a relatively simpler manner and the variational foundation of the phase-field methodprovides provable energy stability and discrete conservation [17]. For these reasons, the phase- field method has attracted more interest in the modeling of the two-phase flow problems in therecent years. The phase-field method considers a diffuse representation of the interface geometry anddescribes the minimization of the free energy functional [18]. The diffuse interface betweenthe two phases is described as a region in which the phases are mixed and store the free energy.The free energy functional E can be written as: E : H (Ω) ∩ L (Ω) → R > , E ( φ ( x , t )) = Z Ω (cid:18) F ( φ ( x , t )) + ε |∇ φ ( x , t ) | (cid:19) d Ω , (1)where Ω is a bounded fluid domain consisting of spatial points x at time t , H (Ω) denotes thespace of square-integrable real-valued functions with square-integrable derivatives on Ω , L (Ω) denotes the function space in which the fourth power of the function is integrable, R > repre-sents the set of non-negative real numbers, φ ( x , t ) is referred to as the order parameter or thephase-field function which indicates the components of the two-phase mixture. The first term F ( φ ( x , t )) in Eq. (1), which is called the bulk or mixing energy, depends on the local compo-sition of the two phases mixture indicated by φ . To invoke the phase separation, a double-well potential is employed: F ( φ ) = ( φ − , in which the minimum bulk energy is attainedwith separated pure phases φ = 1 and φ = − . The second term, which is called interfacial orgradient energy, depends on the composition of the immediate environment indicated by ∇ φ .The gradient energy dictates the interaction between the two phases. The ratio between thebulk energy and the interfacial energy is determined by the interface thickness parameter ε . In the free energy minimization, the bulk energy minimization prefers pure components and sep-arated phases, while the interfacial energy minimization prefers a mixed uniform phase. Using3he interplay between these two effects, the interface thickness parameter ε controls the diffuseinterface geometry.The phase-field methods are generally based on the Cahn-Hilliard and the Allen-Cahnphase-field equations [19, 20], which considers the gradient flow minimizing the free energyfunctional as the cause of the phase-field function evolution. The Cahn-Hilliard equation sat-isfies the mass conservation naturally [21]. However, the equation is a fourth-order partialdifferential equation (PDE), which is cumbersome during the numerical discretization. In con-trast, the Allen-Cahn equation is a second-order convection-diffusion-reaction PDE which hasattractive numerical properties from the implementation standpoint. Although the originalAllen-Cahn equation is not mass-conservative, the conservation property can be realized byadding a Lagrange multiplier [22, 23] or employing an anti-curvature term [24]. The formeris more stable, while the latter is more accurate [25]. Concerning the computational efficiencyand the stability, we employ the Allen-Cahn phase-field equation with a Lagrange multiplierfor solving two-phase flow problems in the current study. In the original Allen-Cahn equation,the evolution of the phase-field function seeks the minimum of the free energy functional: ∂φ∂t = − γ (cid:18) δ E ( φ ) δφ (cid:19) , (2)where γ is the mobility coefficient and δ E ( φ ) δφ = ( F ′ ( φ ) − ε ∇ φ ) represents the variationalderivative of the free energy functional. Eq. (2) can be formulated as the gradient flow of thefree energy functional in L space [26]: ∂∂t E ( φ ) = − γ R Ω (cid:12)(cid:12)(cid:12) δ E ( φ ) δφ (cid:12)(cid:12)(cid:12) d Ω . The mobility coefficientdetermines the intensity of the gradient flow of the free energy functional and controls the speedat which the interface geometry relaxes to the equilibrium profile and shape with minimum freeenergy. For a planar interface in equilibrium, which can be considered as a one-dimensionalcase, the equilibrium interface profile can be solved as: φ eq ( n ) = tanh (cid:16) n √ ε (cid:17) , where n is thecoordinate normal to the interface. The equilibrium interface profile is shown in Fig. 1 (a).The thickness of the diffuse interface is O ( ε ) . When the phase-field method is used for thefluid-fluid interface evolution in the two-phase flow problems, the convection of the flow fieldand the volume conservation must be considered, which leads to the convective form of theconservative Allen-Cahn phase-field equation: ∂φ∂t + u · ∇ φ = − γ (cid:18) δ E ( φ ) δφ − β ( t ) p F ( φ ) (cid:19) , (3)where u represents the convective velocity of the fluid flow, β ( t ) is the time-dependent part of the Lagrange multiplier for mass conservation [23], which is given by β ( t ) = R Ω F ′ ( φ ) d Ω R Ω √ F ( φ ) d Ω . In the convective form of the conservative Allen-Cahn equation, the evolution of the orderparameter is driven by both the convection and the free energy minimization. As the inter-face profile evolves, instead of the equilibrium interface profile, an actual interface profile is formed as a consequence of the interplay between the convective distortion and the free energyminimization. When the finite thickness interface region is subjected to a positive or nega-tive velocity gradient in its normal direction, which represents an extensional or compressionalvelocity field, the interface will be extended or compressed and deviate from the equilibrium4 φ φ eq φ a ConvectivedistortionFree energyminimization ∼ O ( ε )Ω ( φ = −
1) Ω ( φ = 1) xφ (a) φ eq = 0 φ = 0 v ( x , t ) u ( x , t ) (b)Figure 1: Illustrations of the interface dynamics of the convective form of the conservative Allen-Cahn equation:(a) one-dimensional equilibrium interface profile φ eq and the actual interface profile φ a subjected to an extensionalvelocity field, and (b) volume-conserved mean curvature flow velocity v ( x , t ) and the convective velocity u ( x , t ) of the interface φ ( x , t ) = 0 . The free energy minimization described by the equation balances the convectivedistortion with φ a in (a), and induces v ( x , t ) in (b). profile. On the other hand, the deviation from the equilibrium profile increases the free energy. The free energy minimization starts to drive the order parameter back to the equilibrium pro-file. Consequently, an actual profile φ a different from the equilibrium interface profile φ eq isformed due to the interplay, as illustrated in Fig. 1 (a). As mentioned earlier, the intensity of thegradient flow minimizing the free energy functional is controlled by γ . If the free energy min-imization dominates the competition over the convective distortion, the actual interface profile will be kept close to the equilibrium profile, thus achieving interface preservation.During the interface evolution, the volume-preserving mean curvature flow induced by thefree energy minimization disturbs the convection according to the fluid flow velocity. The freeenergy of the interface is closely related to the perimeter or area of the interface [27]. Under theconstraint of the volume conservation, which is the equivalence of the mass conservation witha constant density of each phase, the free energy minimization will drive the actual interfacecontour φ a = 0 towards the contour with minimum perimeter or area. Consequently, theinterface contour convected by the fluid flow velocity relaxes towards a circular or sphericalshape in two or three dimensions, which are denoted as φ eq = 0 [22, 23, 27]. The flow inducedby this process is referred to as the volume-preserving mean curvature flow and illustrated inFig. 1 (b). The velocity of the flow can be derived from the asymptotic analysis [22, 23, 28].For the current formulation, consider the interface Γ φI ( t ) = { x ∈ Ω | φ ( x , t ) = 0 } , the velocityof the volume-preserving mean curvature flow at the interface is given by: v ( x , t ) = γε κ ( x , t ) − | Γ φI ( t ) | Z Γ φI κ ( x , t ) ds ! n φL ( x , t ) , x ∈ Γ φI ( t ) , (4)where v ( x , t ) is the velocity of the volume-preserving mean curvature flow, κ ( x , t ) is definedas κ ( x , t ) = n sd − P i =1 κ i ( x , t ) , n sd being the number of dimensions and κ i ( x , t ) being the prin-ciple curvatures of the interface, | Γ φI ( t ) | is the perimeter or area of the interface, n φL ( x , t ) = ∇ φ/ |∇ φ | is the unit normal vector of the level sets of φ , κ ( x , t ) is defined as positive when n φL ( x , t ) is pointing to the concave side of the interface.5 .3. Related work and contributions The phase-field parameters, namely the mobility coefficient γ and the interface thickness ε , play important roles in the interface-preserving capability and the volume-preserving meancurvature flow velocity. These parameters should be judiciously select to produce a physically- consistent interface behavior. Jacqmin [11] suggested that the gradient flow minimizing thefree energy should properly oppose the convective distortion, while the gradient flow shouldconverge to zero (i.e., the phase-field equation converges to a pure convection equation) as thediffuse interface converges to the sharp interface. According to the order of magnitude analy-sis, a mobility coefficient varying between O ( ε ) and O ( ε ) in the Cahn-Hilliard equation was found to be appropriate. In [29], the magnitude of the total free energy is adjusted dynamicallyto ensure its consistency with the surface tension force coefficient viewed as free energy den-sity. The explicit calculation of the interface length or area was required during the adjustment.In [24], the phase-field propagation equation was formulated for tracking sharp interfaces. Themobility coefficient, the interface thickness parameter, and the maximum convective velocity are combined to construct a non-dimensionalized parameter via a standard explicit finite differ-ence discretization. The parameter was considered purely as a numerical parameter, the impactof which on the interface profile and evolution was studied for stationary and evolving inter-faces. It was found that the increase in the mobility coefficient results in better enforcement ofthe hyperbolic tangent phase-field profile and helps to suppress the instabilities at corners. The mobility coefficient was controlled one order of magnitude below its upper limit given by theCourant-Friedrichs-Levy (CFL) condition to avoid significant discretization errors. In [30], anadditional free energy functional punishing the deviation from the hyperbolic tangent profile isdesigned, which provides a correction term in the Cahn-Hilliard equation via the minimizationprocess. The profile correction term enforces the hyperbolic tangent profile, thereby reducing the interface shrinkage effect, the convective distortion, and improving the surface tension forcecalculation. This approach was further improved in the profile-flux correction [31] and appliedin the turbulent multi-phase flow problems [32].In the current study, we propose an interface-preserving and conservative phase-field methodfor incompressible two-phase flows with a particular emphasis on the accurate surface tension dynamics. A continuum formulation and a systematic approach for determining the phase-fieldparameters γ and ε are presented. The parameters are formulated by directly considering the as-sociated terms in the convective form of the Allen-Cahn equation in a non-dimensional movingorthogonal curvilinear coordinate system. The term representing the effect of the convectivedistortion is identified wherein the convective distortion parameter quantifies the ratio between the convective distortion and the free energy minimization. An interface-preserving conditionfor the parameter enforcing the free energy minimization dominance over the convective dis-tortion is derived. To fulfill the condition, we propose a time-dependent mobility model forcontrolling the RMS convective distortion parameter in the diffuse interface region. Direct re-lationships between the RMS convective distortion parameter and relative interface thickness and surface tension force errors are assessed by numerical simulations of the interface convec-tion problems. By establishing a suitable range of mobility coefficient, the excessive gradientflow minimizing the free energy functional and the resulting spurious volume-preserving meancurvature flow are avoided.The present study builds upon our previous conservative and energy stable variational scheme for the Allen-Cahn and Navier-Stokes system proposed in [17]. The scheme was inte-grated with a mesh adaptivity process in [33], and has been proven to be accurate and stable6or a wide range of fluid-structure interaction problems in the inertia dominate regime withhigh density ratio [34]. In the current work, we further improve the accuracy of the schemein the capillary dominated regime by considering the interface-preserving Allen-Cahn based phase-field model. We employ the model together with the CSF model, where the surfacetension force is transformed into a volume force spread over a few layers of elements. Wediscretize the incompressible Navier-Stokes and Allen-Cahn equations with the finite elementmethod in a fully implicit manner. We maintain the bounded and stable solution of the Allen-Cahn system via the positivity preserving variational (PPV) technique [35] and the coupling between the Allen-Cahn and the Navier-Stokes systems retains second-order accuracy in timedomain [17, 34]. With the aid of a generic 1D bistable diffusion-reaction system in a stretch-ing flow, we first carry out a systematic convergence and verification study of our 1D steadyAllen-Cahn solver based on the PPV technique and the implicit discretization. To demonstratethe interface-preserving formulation, we employ the Allen-Cahn solver for the convection of diffuse interfaces in prescribed incompressible velocity fields for planar and curved situations.We examine the proposed formulation in the two- and three- dimensional rising bubble bench-mark cases through a systematic convergence study. We compare accuracy and convergencewith the sharp interface formulation. Our results show that only when the interface-preservingcapability is improved and the volume-preserving mean curvature flow is decreased simultane- ously, the simulation results will converge to the accurate solution. This requires the reductionof the RMS convective distortion parameter and the interface thickness parameter at the sametime. Finally, we simulate two rising bubbles merging with a free surface with an unstructuredmesh to demonstrate the applicability of the proposed model in practical problems, which hascomplex topological changes of the interface and complex dynamics including bubble-bubble and bubble-free surface interaction.The organization of this paper is as follows: Section 2 presents a mathematical analysisof the diffuse interface profile, wherein the convective distortion parameter and the interface-preserving condition are identified. The time-dependent mobility model is proposed accordingto the interface-preserving condition. Section 3 describes the implementation of the varia- tional formulation for the interface-preserving conservative Allen-Cahn-Navier-Stokes systemwith the time-dependent mobility model. Section 4 verifies our implementation through a 1Dbistable convection-diffusion-reaction system in a stretching flow and provides a numerical as-sessment of the errors associated with the convective distortion parameter for a planar and acurved interface. Two- and three-dimensional rising bubble cases are investigated in Section
2. Interface-preserving phase-field formulation
In this section, we present the continuum formulation of the time-dependent mobility modelfor preserving the hyperbolic tangent profile. The convective form of the Allen-Cahn equationis directly analyzed in a non-dimensional moving orthogonal curvilinear coordinate system.The term representing the influence of the convective distortion is identified in the governingequation. The magnitude of the term depends on a non-dimensional parameter, which we refer to as the convective distortion parameter. An interface-preserving condition is derived for the7onvective distortion parameter to preserve the interface profile. The time-dependent mobilitymodel is proposed based on the interface-preserving condition.
Following the original work of [20], we describe the evolution of the two-phase interface inan orthogonal curvilinear coordinate system. This allows us to simplify the governing equationutilizing the property that the level sets of the order parameter are parallel to the interface. Con-sider a physical domain Ω × ]0 , T [ with spatial coordinates x and temporal coordinate t . Theboundary of the computational domain Γ is decomposed as Γ = Γ φD ∪ Γ φH , where Γ φD and Γ φH denote the Dirichlet and Neumann boundaries for the order parameter respectively. A diffuseinterface that separates the immiscible two-phase fluids defined on Ω is indicated by the orderparameter φ ( x , t ) . The diffuse interface is convected by a velocity field u ( x , t ) . In the orthog-onal curvilinear coordinate system, the spatial coordinates are given by x = ( n, τ , τ ) , where n is the coordinate of the axis normal to the level sets of φ , and the rest two coordinates τ and τ are the coordinates of the axes which are tangential to the level sets of φ . In this coordinatesystem, the convection of a diffuse interface is given by the following initial boundary valueproblem based on the Allen-Cahn equation: ∂φ∂t + u · ∇ φ = − γ (cid:0) F ′ ( φ ) − ε ∇ φ (cid:1) , on Ω ,φ = φ D , ∀ x ∈ Γ φD , n φ Γ · ∇ φ = 0 , ∀ x ∈ Γ φH ,φ (cid:12)(cid:12) t =0 = φ , on Ω , (5)where F ′ ( φ ) = φ − φ is the derivative of the double-well potential with respect to φ , n φ Γ repre-sents the unit vector normal to the boundary of the computational domain and φ represents theinitial condition for the order parameter. The velocity in the orthogonal curvilinear coordinatesystem can be expressed as: u = u n n φL + u τ τ + u τ τ , (6)where n φL , τ and τ are the unit vectors in the normal and two tangential directions of the level sets of φ respectively, and u n , u τ , u τ are the corresponding velocity components.Assume that the normal profile of the interface is almost the same everywhere on the interface.Therefore the derivatives of the order parameter in the tangential directions are negligible. Withthis assumption, the spatial derivatives of the order parameter can be calculated as: ∇ φ = ∂φ∂n n φL , ∇ φ = ∇ · ∇ φ = ∂ φ∂n + ∂φ∂n ∇ · n φL . (7)Notice that ∇ · n φL = − κ , where κ is the summation of principle curvatures of the interface.Substituting Eqs. (6) and (7) into the first equation of Eq. (5), the convective form of the Allen-Cahn equation in the orthogonal curvilinear coordinate system is given by: ∂φ∂t + u n ∂φ∂n = − γ (cid:18) F ′ ( φ ) − ε (cid:18) ∂ φ∂n − κ ∂φ∂n (cid:19)(cid:19) . (8)8 ( x , φ ( x , t ) u ( x , t ) ˜ n m ˜ τ ˜ τ ∼ O ( ε ) ∼ O (1) Figure 2: Schematic diagram of the diffuse interface in non-dimensional moving orthogonal curvilinear coordinatesystem. The coordinate system ( ˜ n m , ˜ τ and ˜ τ ) is attached on the moving interface indicated by φ ( x , t ) convectedin the velocity field u ( x , t ) . The thickness of the diffuse interface in the coordinate system is ∼ O (1) due to thenon-dimensionalization. Now we write the convective Allen-Cahn equation in a non-dimensional moving orthogonalcurvilinear coordinate system. We non-dimensionalize the coordinate system using the inter-face thickness parameter. For the convenience of analyzing the interface distortion due to theconvective velocity difference in the diffuse interface region, we translate the coordinate sys-tem with the interface. As a result, the relative convective velocity to the interface, which leadsto the convective distortion, appears explicitly in the governing equation. To begin with, wecarry out the non-dimensionalization of the coordinate system by denoting the dimensionlesscoordinates as ˜ n, ˜ τ , ˜ τ , which are non-dimensionalized with respect to ε : ˜ n = n/ε, ˜ τ = τ /ε, ˜ τ = τ /ε. (9)Non-dimensionalizing the spatial coordinates in Eq. (8) accordingly, we obtain: ∂φ∂t + u n ∂φ∂ ˜ n ε = − γ (cid:18) F ′ ( φ ) − (cid:18) ∂ φ∂ ˜ n − εκ ∂φ∂ ˜ n (cid:19)(cid:19) . (10)Assume that the principal radii of the interface are large compared to the interface thickness,which leads to κ ≪ /ε . With this assumption, the last term in Eq. (10) can be neglected: ∂φ∂t + u n ∂φ∂ ˜ n ε = − γ (cid:18) F ′ ( φ ) − ∂ φ∂ ˜ n (cid:19) . (11)This completes the non-dimensionalization.Next, we write Eq. (11) in a coordinate system which translates with the interface. We definethe coordinates of the interface as Γ φI ( t ) = { ( n ( t ) , τ , τ ) ∈ Ω | φ ( n , τ , τ , t ) = 0 } . Thecoordinate transformation to the moving coordinate system can be written as: n m = n − n ( t ) , ˜ n m = ˜ n − ˜ n ( t ) , ˜ n ( t ) = n ( t ) /ε, (12)9here n m represents the normal coordinate in the moving coordinate system, ˜ n m is the non-dimensionalized n m with respect to ε , ˜ n ( t ) is the non-dimensional normal coordinate of theinterface. The coordinate system is illustrated in Fig. 2. In the moving coordinate system, the temporal and spatial derivatives of the order parameterbecome: ∂φ (˜ n m , ˜ τ , ˜ τ , t ) ∂t = ∂φ∂t − ∂φ∂ ˜ n m d ˜ n ( t ) dt , ∂φ∂n = ∂φ∂n m , ∂φ∂ ˜ n = ∂φ∂ ˜ n m , ∂ φ∂ ˜ n = ∂ φ∂ ˜ n m . (13)Replacing the temporal and spatial derivatives in Eq. (11) with Eq. 13, we have: ∂φ∂t − ∂φ∂ ˜ n m d ˜ n ( t ) dt + u n (˜ n m , ˜ τ , ˜ τ , t ) ε ∂φ∂ ˜ n m = − γ (cid:18) F ′ ( φ ) − ∂ φ∂ ˜ n m (cid:19) . (14)Notice that ˜ n ( t ) is the non-dimensional normal coordinate of the interface, the time derivativeof which in Eq. (14) gives the normal velocity of the interface in the non-dimensional coor-dinate system. To get the normal velocity at the interface, we substitute φ = 0 into Eq. (11): ∂φ∂t (cid:12)(cid:12)(cid:12)(cid:12) φ =0 + u n (˜ n ( t ) , ˜ τ , ˜ τ , t ) ∂φ∂ ˜ n (cid:12)(cid:12)(cid:12)(cid:12) φ =0 ε = − γ − ∂ φ∂ ˜ n (cid:12)(cid:12)(cid:12)(cid:12) φ =0 ! . (15)For the hyperbolic tangent profile, ∂ φ∂ ˜ n = 0 at φ = 0 . Assuming that this is approximatelysatisfied when the convective distortion is not significant, thus the right-hand side of Eq. (15)is negligible: ∂φ∂t (cid:12)(cid:12)(cid:12)(cid:12) φ =0 + u n (˜ n ( t ) , ˜ τ , ˜ τ , t ) ε ∂φ∂ ˜ n (cid:12)(cid:12)(cid:12)(cid:12) φ =0 = 0 . (16)The velocity of the interface can be identified from the convection equation ( ) as: d ˜ n ( t ) dt = u n (˜ n ( t ) , ˜ τ , ˜ τ , t ) ε . (17)Substituting Eq. (17) into Eq. (14), rewriting the velocity in the moving coordinate systemand non-dimensionalize the equation with respect to γ , we get the governing equation for theinterface profile: ∂φ∂ ˜ t + (cid:18) u n (˜ n m , ˜ τ , ˜ τ , ˜ t ) − u n (0 , ˜ τ , ˜ τ , ˜ t ) γε (cid:19) ∂φ∂ ˜ n m = − (cid:18) F ′ ( φ ) − ∂ φ∂ ˜ n m (cid:19) , (18)where ˜ t = tγ is the non-dimensional time. As the second term in Eq. (18) goes to zero, theequation recovers to the Allen-Cahn equation, which gives the well-known hyperbolic tangentprofile. The assumptions used in the derivation that the normal interface profile is almostidentical on the interface and ∂ φ∂ ˜ n = 0 is approximately satisfied at φ = 0 are valid. As mentioned earlier, a non-zero second term of Eq. (18) causes the deviation from thehyperbolic tangent profile due to convection. We refer to the term as the convective distortionterm. Multiplying and dividing the term by ˜ n m , Eq. (18) becomes: ∂φ∂ ˜ t + (cid:18) u n (˜ n m , ˜ τ , ˜ τ , ˜ t ) − u n (0 , ˜ τ , ˜ τ , ˜ t )˜ n m (cid:19) γε ˜ n m ∂φ∂ ˜ n m = − (cid:18) F ′ ( φ ) − ∂ φ∂ ˜ n m (cid:19) . (19)10uppose that u n is continuous on [0 , ˜ n m ] and differentiable on (0 , ˜ n m ) when < ˜ n m , or u n iscontinuous on [˜ n m , , and differentiable on (˜ n m , when ˜ n m < , according to the mean valuetheorem [36], there exists ˜ n s ∈ (0 , ˜ n m ) or ˜ n s ∈ (˜ n m , respectively such that: (cid:18) u n (˜ n m , ˜ τ , ˜ τ , ˜ t ) − u n (0 , ˜ τ , ˜ τ , ˜ t )˜ n m (cid:19) = ∂u n ∂ ˜ n m (˜ n s , ˜ τ , ˜ τ , ˜ t ) . (20)With Eq. (20), Eq. (19) can be rewritten as: ∂φ∂ ˜ t + (cid:18) ∂u n ∂ ˜ n m (˜ n s , ˜ τ , ˜ τ , ˜ t ) (cid:19) γε ˜ n m ∂φ∂ ˜ n m = − (cid:18) F ′ ( φ ) − ∂ φ∂ ˜ n m (cid:19) , (21)where ˜ n s ∈ (0 , ˜ n m ) when < ˜ n m and ˜ n s ∈ (˜ n m , when ˜ n m < . Transform the par-tial derivative of the normal velocity back to the dimensional spatial coordinate system usingEqs. (9) and (13): ∂φ∂ ˜ t + (cid:18) ∂u n ∂n ( n s , τ , τ , ˜ t ) (cid:19) γ ˜ n m ∂φ∂ ˜ n m = − (cid:18) F ′ ( φ ) − ∂ φ∂ ˜ n m (cid:19) , (22)where n s ∈ (0 , n m ) when < n m and n s ∈ ( n m , when n m < .For the convenience of notation, we replace the notation of the non-dimensional time ˜ t with t and denote the normal velocity gradient in the normal direction as: ζ ( x , t ) = ∂u n ∂n ( x , t ) , (23)which can be considered as the intensity of the convective distortion (see Appendix A for de-tailed explanation). Remark 1.
The magnitude of ζ ( x , t ) usually increases with the increase in the principal cur- vatures of the interface. Because high principal curvatures lead to a large surface tension force.When the force is distributed to the diffuse interface region by a Dirac delta function, the vi-olent variation of the surface tension force term across the diffuse interface gives rise to highvelocity gradients, and furthermore causes the large magnitude of ζ ( x , t ) . We next define a non-dimensional parameter, which we refer to as the convective distortionparameter: ξ ( x , t ) = ζ ( x , t ) (cid:14) γ. (24)With the notations in Eqs. (23) and (24), the governing equation for the interface profile Eq. (22)becomes: ∂φ∂t + ξ ( n s , τ , τ , t )˜ n m ∂φ∂ ˜ n m = − (cid:18) F ′ ( φ ) − ∂ φ∂ ˜ n m (cid:19) , (25)where n s ∈ (0 , n m ) when < n m and n s ∈ ( n m , when n m < .As the convective distortion term goes to zero, the interface profile approaches the hyper-bolic tangent profile. Thus, the interface-preserving condition is given by: (cid:12)(cid:12)(cid:12)(cid:12) ξ ( n s , τ , τ , t )˜ n m ∂φ∂ ˜ n m (cid:12)(cid:12)(cid:12)(cid:12) η , (26)11here η is a desired upper bound for the magnitude of the convective distortion term. Applythe CauchySchwarz inequality for the convective distortion term in the condition (26) as: (cid:12)(cid:12)(cid:12)(cid:12) ξ ( n s , τ , τ , t )˜ n m ∂φ∂ ˜ n m (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) ξ ( n s , τ , τ , t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ n m ∂φ∂ ˜ n m (cid:12)(cid:12)(cid:12) η , (27)where n s ∈ (0 , n m ) when n m > and n s ∈ ( n m , when n m < . For the hyperbolic tangentprofile φ (˜ n m ) = tanh(˜ n m / √ , | ˜ n m ∂φ∂ ˜ n m | approaches zero outside the diffuse interface region,which satisfies the condition (27). Inside the diffuse interface region, | ˜ n m ∂φ∂ ˜ n m | is bounded by aconstant. Denoting the constant as η , we have: (cid:12)(cid:12)(cid:12)(cid:12) ˜ n m ∂φ∂ ˜ n m (cid:12)(cid:12)(cid:12)(cid:12) η , x ∈ Γ φDI ( t ) , (28)where Γ φDI ( t ) represents the diffuse interface region at time t . We assume that inequality (28)is still valid when the interface profile is close to the hyperbolic tangent profile. Thus, theinterface-preserving condition becomes: | ξ ( x , t ) | η, x ∈ Γ φDI ( t ) (29)where η = η /η is a desired upper bound for the convective distortion parameter. Since η is aconstant once the diffuse interface region is defined, η → as η → . In other words, the mag-nitude of the convective distortion term solely depends on ξ and decreases with the reductionof the absolute value of the convective distortion parameter, which leads to the convergence of the interface profile to the hyperbolic tangent profile. Remark 2.
The diffuse interface region defines the spatial domain where the convective dis-tortion is considered. A large diffuse interface region includes the distortion away from theinterface, which may lead to an overestimation. The overestimation provides better interface- preserving capability but may induce a relatively larger volume-preserving mean curvatureflow. In contrast, a small diffuse interface region may underestimate the interface distortiongiving rise to a smaller volume-preserving mean curvature flow but weakening the interface-preserving capability. In the current study, we define the diffuse interface region as the regionwhere of the variation of φ occurs: Γ φDI ( t ) = { ( x , t ) || φ ( x , t ) | . } . For the hyperbolic tangent profile, the thickness of the region in the normal direction is . ε . To satisfy the interface-preserving condition, one needs to adjust the mobility dynamicallyaccording to the normal velocity gradient in the normal direction. The interface-preservingcondition (29) can be written as: γ > | ζ ( x , t ) | η , x ∈ Γ φDI ( t ) . (30)We notice that the right-hand side of inequality (30) varies in both space and time. To maintainthe inequality, it is natural to consider a mobility model with spatial and temporal dependence.However, a mobility model that prohibits the variation in the normal direction while allows12ariations in tangential directions poses challenges in its construction. Consequently, we con-sider a time-dependent mobility model in the current study. The mobility coefficient is takenas a constant throughout the computational domain at each time instance, while it is allowed tochange as time evolves. This requires a projection at time t from the spatially varying | ζ ( x , t ) | to a real-valued γ ( t ) : F : Ω → R > , γ ( t ) = 1 η F ( | ζ ( x , t ) | ) . (31)In the current study, we employ the RMS value in the diffuse interface region for the projection,which relaxes the condition (30) in an average sense: F ( ϕ ( x , t )) = s R ( ϕ ( x , t )) d Ω R d Ω , x ∈ Γ φDI ( t ) . (32)We refer to η as the RMS convective distortion parameter in the rest of the paper. As de-rived in Appendix B, the frame independent form of the time-dependent mobility model can beexpressed as: γ ( t ) = 1 η F (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ∇ φ · ∇ u · ∇ φ |∇ φ | (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) . (33)The frame independent form facilitates its numerical implementation in Cartesian coordinatesystem introduced in the next section. Remark 3.
By taking the RMS function to perform the projection, the condition (30) is relaxedin an average sense. Because the mobility is calculated according to the RMS value of | ζ ( x , t ) | ,at the location where | ζ ( x , t ) | exceeds the RMS value, the condition (30) is violated. As men-tioned in Remark 1, the violation usually happens in the region with high principle curvaturesand singularities of the interface. This allows the merging and breaking-up of the interface, where keeping the hyperbolic tangent profile is no longer required. Other projection methodscan be employed for different considerations and requirements.
3. Variational interface-preserving conservative Allen-Cahn-Navier-Stokes formulation
In this section, we present an variational implementation of the interface-preserving con- servative phase-field formulation. For the sake of completeness, we describe the governingequations of the two-phase flow modeling, viz., the incompressible Navier-Stokes equationsand the conservative Allen-Cahn equation with the proposed time-dependent mobility model.We begin with the strong form of the equations and then project them into finite element spaceas the semi-discrete variational form. Specifically, we describe the discretization of the time- dependent mobility model. The section is closed with the coupled linearized matrix form ofthe variational discretization.
Consider a domain Ω × ]0 , T [ consisting of the spatial points x at time t . The boundary ofthe domain, Γ can be decomposed in two ways, Γ = Γ f D ∪ Γ f H and Γ = Γ φD ∪ Γ φH , where Γ f D and Γ f H denote the Dirichlet and Neumann boundaries for the Navier-Stokes equations respectively,13hile Γ φD and Γ φH denote the same for the Allen-Cahn counterpart respectively. The diffuseinterface region between the two-phases is denoted as Γ φDI ( t ) . The one-fluid formulation forthe two-phase incompressible and immiscible fluids system with the boundary conditions isgiven as: ρ ∂ u ∂t + ρ u · ∇ u = ∇ · σ + sf + b , on Ω , ∇ · u = 0 , on Ω , u = u D , ∀ x ∈ Γ f D , σ · n f = h , ∀ x ∈ Γ f H , u = u , on Ω(0) , (34) ∂φ∂t + u · ∇ φ = − γ ( t ) (cid:0) F ′ ( φ ) − ε ∇ φ − β ( t ) p F ( φ ) (cid:1) , on Ω ,γ ( t ) = 1 η F ( | ζ ( x , t ) | ) x ∈ Γ φDI ( t ) φ = φ D , ∀ x ∈ Γ φD , ∇ φ · n φ Γ = 0 , ∀ x ∈ Γ φH ,φ (cid:12)(cid:12) t =0 = φ , on Ω(0) , (35)where Eq. (34) and Eq. (35) represent the Navier-Stokes and Allen-Cahn equations respectively.In the Navier-Stokes equations, ρ is the density of the fluid, u represents the fluid velocitydefined for each spatial point x in Ω , b is the body force on the fluid such as gravity ( b = ρ g ) , g being the acceleration due to gravity, u D and h denote the boundary conditions at the Dirichletand Neumann boundaries respectively, n f is the unit outward normal to the Neumann boundaryand u represents the initial velocity field at t = 0 . The Cauchy stress tensor for a Newtonianfluid is given as: σ = − p I + T , T = 2 µ ǫ ( u ) , ǫ ( u ) = 12 (cid:2) ∇ u + ( ∇ u ) T (cid:3) , (36)where p is the pressure field, T and ǫ represent the shear stress tensor and the fluid strain ratetensor respectively and µ denotes the dynamic viscosity of the fluid. The physical parametersof the fluid such as ρ and µ vary with the evolution of the interface indicated by order parameter φ : ρ ( φ ) = 1 + φ ρ + 1 − φ ρ , (37) µ ( φ ) = 1 + φ µ + 1 − φ µ , (38)where ρ i and µ i are the density and dynamic viscosity of the i th phase of the fluid respectively.The surface tension force sf is modeled by the continuum surface force (CSF) model [9], inwhich it is reformulated as a volumetric source term with a Dirac Delta function utilizing thegradient of the phase indicator φ . Several forms of sf ( φ ) have been used in the literature which14re reviewed in [37, 38]. In this study, we employ the following definition [39]: sf ( φ ) = σκ n φL δ S = σ ∇ · (cid:16)(cid:16) I − n φL ⊗ n φL (cid:17) δ S (cid:17) = σα sf ε ∇ · (cid:0) |∇ φ | I − ∇ φ ⊗ ∇ φ (cid:1) (39)where σ is the surface tension coefficient, δ S = εα sf |∇ φ | is the Dirac delta function at theinterface, α sf = 3 √ / is a constant derived by the property of the Dirac delta function, κ being the summation of the principle curvatures of the interface and n φL = ∇ φ/ |∇ φ | denotesthe normal vector of the level sets of φ .On the other hand, in Eq. (35), ε is the interface thickness parameter, γ ( t ) is the time-dependent mobility, F ( φ ) is the double-well potential, η is the RMS convective distortionparameter and n φL is the unit normal vector of the level sets of the order parameter φ . Thevalue of the order parameter at the Dirichlet boundary is denoted by φ D , the initial condition isrepresented by φ and n φ Γ denotes the unit outward normal to the Neumann boundary where azero flux condition is satisfied. The mass conservation is enforced in the Allen-Cahn equationby a Lagrange multiplier β ( t ) p F ( φ ) where β ( t ) = R Ω F ′ ( φ )dΩ / R Ω p F ( φ )dΩ , F ′ ( φ ) is thederivative of the energy potential with respect to the order parameter. The Allen-Cahn equationcan be transformed into a convection-diffusion-reaction equation as follows: ∂ t φ + u · ∇ φ − γ ( t )(ˆ k ∇ φ − ˆ sφ + ˆ f ) = 0 on Ω f , (40)where u , ˆ k , ˆ s and ˆ f are the convective velocity, modified diffusion coefficient, modified reac-tion coefficient and the modified source respectively which are defined in [17]. In this subsection, we present the semi-discrete variational form of the Navier-Stokes-Allen-Cahn (NS-AC) system, which has been described earlier. We employ the generalized- α technique [40] for the temporal discretization which enables a user-controlled high frequencydamping desirable for coarse discretizations in space and time. The following expressions areemployed for the temporal discretization of the Navier-Stokes equations: u n+1 = u n + ∆ t∂ t u n + ς ∆ t ( ∂ t u n+1 − ∂ t u n ) , (41) ∂ t u n+ α m = ∂ t u n + α m ( ∂ t u n+1 − ∂ t u n ) , (42) u n+ α = u n + α ( u n+1 − u n ) , (43)where α , α m and ς are the generalized- α parameters which are dependent on the user-definedspectral radius ρ ∞ . The time step size is denoted by ∆ t and ∂ t denotes the partial differentiationwith respect to time. Similar expressions can be written for the Allen-Cahn equation as well.Suppose S h u , S h p and S h φ denote the space of trial solution such that: S h u = (cid:8) u h | u h ∈ ( H (Ω)) d , u h = u D on Γ f D (cid:9) , (44) S h p = (cid:8) p h | p h ∈ L (Ω) (cid:9) , (45) S h φ = (cid:8) φ h | φ h ∈ H (Ω) , φ h = φ D on Γ φD (cid:9) , (46)15here ( H (Ω)) d denotes the space of square-integrable R d -valued functions with square-integrablederivatives on Ω and L (Ω) is the space of the scalar-valued functions that are square-integrableon Ω . Similarly, we define V h ψ , V h q and V h φ as the space of test functions such that: V h ψ = (cid:8) ψ h | ψ h ∈ ( H (Ω)) d , ψ h = on Γ f D (cid:9) , (47) V h q = (cid:8) q h | q h ∈ L (Ω) (cid:9) , (48) V h φ = (cid:8) ˆ w h | ˆ w h ∈ H (Ω) , ˆ w h = 0 on Γ φD (cid:9) . (49)The variational statement of the combined NS-AC system can be written as:find [ u h ( t n+ α ) , p h ( t n+1 ) , φ h ( t n+ α )] ∈ S h u × S h p × S h φ such that ∀ [ ψ h , q h , ˆ w h ] ∈ V h ψ × V h q × V h φ for the incompressible NS equations Z Ω ρ ( φ )( ∂ t u h + u h · ∇ u h ) · ψ h dΩ + Z Ω σ h : ∇ ψ h dΩ + Z Ω α sf σε (cid:0) |∇ φ h | I − ∇ φ h ⊗ ∇ φ h (cid:1) : ∇ ψ h dΩ | {z } Surface tension force + n el X e=1 Z Ω e τ m ρ ( φ ) ( ρ ( φ ) u h · ∇ ψ h + ∇ q h ) · R m dΩ e + Z Ω q h ( ∇ · u h )dΩ + n el X e=1 Z Ω e ∇ · ψ h τ c ρ ( φ ) R c dΩ e − n el X e=1 Z Ω e τ m ψ h · ( R m · ∇ u h )dΩ e − n el X e=1 Z Ω e ∇ ψ h ρ ( φ ) : ( τ m R m ⊗ τ m R m )dΩ e = Z Ω b ( t n+ α ) · ψ h dΩ + Z Γ H h · ψ h dΓ , (50)and for the Allen-Cahn equation: Z Ω (cid:18) ˆ w h ∂ t φ h + ˆ w h (cid:0) u h · ∇ φ h (cid:1) + γ (cid:0) t n + α (cid:1)| {z } Dynamic mobility (cid:0) ∇ ˆ w h · (ˆ k ∇ φ h ) + ˆ w h ˆ sφ h − ˆ w h ˆ f (cid:1)(cid:19) dΩ+ n el X e=1 Z Ω e (cid:18)(cid:16) u h · ∇ ˆ w h (cid:17) τ φ (cid:16) ∂ t φ h + u h · ∇ φ h − γ (cid:0) t n + α (cid:1)| {z } (cid:0) ∇ · (ˆ k ∇ φ h ) − ˆ sφ h + ˆ f (cid:1)(cid:17)(cid:19) dΩ e + n el X e=1 Z Ω e χ |R ( φ h ) ||∇ φ h | k add s ∇ ˆ w h · (cid:18) u h ⊗ u h | u h | (cid:19) · ∇ φ h dΩ e + n el X e=1 Z Ω e χ |R ( φ h ) ||∇ φ h | k add c ∇ ˆ w h · (cid:18) I − u h ⊗ u h | u h | (cid:19) · ∇ φ h dΩ e = 0 , (51)where the terms with under brackets representing the fluid-fluid interface dynamics are centralto the current study, R m , R c and R ( φ h ) denote the element-wise residuals for the momentum, continuity and the Allen-Cahn equations, respectively.In Eq. (50), the terms in the first line represent the Galerkin projection of the momentumequation in the test function space ψ h and the second line comprises of the Petrov-Galerkinstabilization term for the momentum equation. The third line denotes the Galerkin projectionand stabilization terms for the continuity equation and the terms in the fourth line are derived16ia approximation of fine scale velocity on the element interiors based on multi-scale argument[41, 42]. The terms in the final line are the Galerkin projection of the body force and Neumannboundary condition. On the other hand, in Eq. (51), the first line is the Galerkin projection ofthe transient, convection, diffusion, reaction and source terms, the second line represents theStreamline-Upwind Petrov-Galerkin stabilization and the third line depicts the PPV terms thatare derived for the multi-dimensional convection-diffusion-reaction equation via satisfaction ofthe positivity condition at the element matrix level [35]. Several test cases have been performedto assess the effectiveness of this PPV technique in [35]. The details of the derivation of theadded diffusions k add s , k add c and χ can be found in [35], which are given for the present contextby [17]: χ = 2 | ˆ s | h + 2 | u h | , (52) k add s = max (cid:26) || u h | − τ φ | u h | ˆ s | h − (ˆ k + τ φ | u h | ) + ˆ sh , (cid:27) , (53) k add c = max (cid:26) | u h | h − ˆ k + ˆ sh , (cid:27) , (54)where | u h | is the magnitude of the convective velocity and h is the characteristic element lengthdefined in [35]. The stabilization parameters τ m , τ c and τ φ in Eqs. (50) and (51) are given by[43, 44]: τ m = (cid:20)(cid:18) t (cid:19) + u h · Gu h + C I (cid:18) µ ( φ ) ρ ( φ ) (cid:19) G : G (cid:21) − / , τ c = 1tr( G ) τ m , (55) τ φ = (cid:20)(cid:18) t (cid:19) + u h · Gu h + 9ˆ k G : G + ˆ s (cid:21) − / . (56)where C I is a constant derived from the element-wise inverse estimates [45], G is the elementcontravariant metric tensor and tr( G ) is the trace of the contravariant metric tensor. This sta-bilization in the variational form circumvents the Babu ˇs ka-Brezzi condition that is required tobe satisfied by any standard mixed Galerkin method [46]. We present the discrete form of the time-dependent mobility model in this subsection. Thestrong form of the time-dependent mobility model is given by Eq. (33): γ ( t ) = 1 η F (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ∇ φ · ∇ u · ∇ φ |∇ φ | (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) , (57)where F ( ϕ ( x , t )) = q R ( ϕ ( x ,t )) d Ω R d Ω , x ∈ Γ φDI ( t ) . In Eq. (57), F (cid:16)(cid:12)(cid:12)(cid:12) ∇ φ ·∇ u ·∇ φ |∇ φ | (cid:12)(cid:12)(cid:12)(cid:17) can be approxi-mated as the RMS of |∇ φ · ∇ u · ∇ φ/ |∇ φ | | at all the nodes located inside the diffuse interfaceregion Γ φDI ( t ) . The nodal value of |∇ φ · ∇ u · ∇ φ/ |∇ φ | | is calculated as follows.The nodal value of u h ( t n+ α ) , φ h ( t n+ α ) is used to interpolate the ∇ u h ( t n+ α ) and ∇ φ h ( t n+ α ) at the quadrature points, and L -projection is used to project the value on the quadrature pointsback to the nodes inside the diffuse interface region [47]. If the node lies outside the diffuse17nterface region, the value is assigned to be zero. For example, for the node p , we have: (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:18) ∇ φ h · ∇ u h · ∇ φ h |∇ φ h | (cid:19) p (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) P e R Ω e N p ( ∇ φ h · ∇ u h · ∇ φ h ) / |∇ φ h | d Ω e P e R Ω e N p d Ω e (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) if | φ p | . , | φ p | > . , where N p represents the shape function at node p . The discrete form of the time-dependentmobility model is given by: γ ( t n+ α ) = 1 η vuut n DI n DI X p =1 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:18) ∇ φ h · ∇ u h · ∇ φ h |∇ φ h | (cid:19) p (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (58)where n DI denotes the number of the nodes lying inside the diffuse interface region. In this subsection, we present the coupled linearized matrix form of the variationally dis-cretized two-phase flow equations. Employing the Newton-Raphson linearization technique,the coupled two-phase fluid system can be expressed in terms of the solution increments forvelocity, pressure and order parameter ( ∆ u , ∆ p and ∆ φ respectively) as: K Ω G Ω D Ω − G T Ω C Ω G AC K AC ∆ u ∆ p ∆ φ = R m R c R ( φ ) (59)where K Ω is the stiffness matrix of the momentum equation consisting of transient, convection,viscous and Petrov-Galerkin stabilization terms, G Ω is the gradient operator, G T Ω is the diver-gence operator for the continuity equation and C Ω is the stabilization term for cross-couplingof pressure terms. On the other hand, D Ω consists of the terms in the momentum equation which depend on the phase-indicator φ , G AC is the velocity coupled term in the Allen-Cahnequation and K AC is the left-hand side stiffness matrix for the Allen-Cahn equation comprisingof transient, convection, diffusion, reaction and positivity preserving stabilization terms. Here, R m , R c and R ( φ ) represent the weighted residuals of the variational forms in Eqs. (50-51).The two-phase flow system in Eq. (59) is decoupled into two subsystems: Navier-Stokesand Allen-Cahn solves, for which the linear system of equations can be summarized as: " K Ω G Ω − G T Ω C Ω ∆ u ∆ p ) = ( R m R c ) (60) (cid:2) K AC (cid:3) (cid:8) ∆ φ (cid:9) = (cid:8) R ( φ ) (cid:9) (61)Note that the cross-coupling terms between the Navier-Stokes and the Allen-Cahn equations ( D Ω and G AC ) are not present in the decoupled form. We solve the decoupled system ina partitioned-block iterative manner which leads to flexibility and ease in its implementationto the existing variational solvers. The linear systems (Eqs. (60) and (61)) are solved by theGeneralized Minimal Residual (GMRES) algorithm proposed by [48]. The algorithm relies on18rylov subspace iteration and modified Gram-Schmidt orthogonalization. Instead of construc- tion of the left-hand side matrices explicitly, we only construct the required matrix-vector prod-ucts of each block matrix in the GMRES solver. Detailed algorithmic steps for the partitionedcoupling of the fully implicit solutions of the conservative Allen-Cahn and the incompressibleNavier-Stokes equations can be found in [17]. The stability and robustness of the partitioneddecoupled system has been demonstrated for a broad range of problems involving high-density and viscosity ratios, high Reynolds number and complex topological changes over unstructuredmeshes [17, 34].
4. Interface convection problem
In this section, we first verify the convergence and accuracy of our fully-implicit finiteelement solver by simulating a bistable steady convection-diffusion-reaction system in a one- dimensional stretching flow. We then turn our attention to the convection of a planar interfaceand a curved interface in prescribed velocity fields. The convective distortion of the diffuseinterface is quantified by the relative interface thickness and surface tension force errors. Thedependence of the errors on the convective distortion parameter ξ is assessed systematically. For simplification, we consider 1D steady-state Allen-Cahn solution of the interface profileEq. (25) for constant convective distortion parameter: ξ ˜ x dφd ˜ x = − ( φ − φ − d φd ˜ x ) , (62)where ˜ x = x/ε is the non-dimensionalized coordinate. The above interface profile equa-tion for a phase-field function can be considered as a special case of a generic bistable steadyconvection-diffusion-reaction system in a stretching flow, which can be written as: − Sx dφdx = D d φdx + Rφ ( φ − A − φ ) , (63)where S is the stretching rate, D is the diffusion coefficient, R is the reaction rate, and A isa parameter determines the unstable equilibrium phase separating the two stable equilibriumphases with minimum bulk energy. In [49], Eq. (63) has been solved semi-analytically in thelarge Damk¨ohler number limit defined as Da = R/S ≫ with a fixed stretching rate S = 1 .The solutions of Eq. (63) include the plateau-like solution and the pulse-like solution. Sincethe plateau-like solution can be considered as a solution composed of two diffuse interfacessubjected to the convective distortion (as shown in Fig. 4 (a)), it is used for the verification ofour implicit finite element Allen-Cahn solver. The plateau-like solution takes the form: φ ( x ) = 12 f [tanh ( w ( x + v )) − tanh ( w ( x − v ))] , where f represents the height of the “plateau”, w is referred to as the inverse width of thediffuse interface which is inversely proportional to the width of the diffuse interface, and v isthe half-width of the “plateau”. They are given semi-analytically by: f ∼ , w ∼ r Da D , v ∼ √ DDa (0 . − A ) . (64)19o form a discrete system which is consistent with Eq. (63), we directly prescribe the ve-locity as u = − Sx without solving the Navier-Stokes equations. The Lagrange parameter is setto be zero and the reaction term is adjusted accordingly. We consider a one-dimensional com-putational domain x ∈ [ − L, L ] . The computational domain is discretized by a uniform meshof grid size h . A zero flux boundary condition is applied on the left and the right boundaries.The initial condition is specified as: φ ( x ) = 12 (cid:18) tanh (cid:18) √ x + 5) (cid:19) − tanh (cid:18) √ x − (cid:19)(cid:19) . (65)The stretching rate and the diffusion coefficient are taken as S = 1 and D = 1 respectively. Thereaction rates R ∈ [20 , corresponding to Da ∈ [20 , are considered for the verificationpurpose. The unstable equilibrium phase is set to be A = 0 . and the time step is taken as ∆ t = 0 . for the marching of steady state solution.We perform a systematic convergence study for the final time when the steady state solution is reached, the length of the computational domain and the grid size. In the convergence study,we set the tolerance to − for both the linear GMRES and the nonlinear Newton solvers.To quantify the error in the convergence study, we define the relative error in L norm as e = || φ − φ ref || / || φ ref || , where || · || denotes the L norm of the vector. We first study thefinal time for the case Da = 140 , h = 0 . and L = 15 . The error level of the steady solution at t = 36 reaches e = 1 . × − while considering the solution at t = 40 as the reference.Hence we consider the solution at t = 40 as the fully-converged steady state solution. We nextinvestigate the domain length L to ensure that the zero flux boundary condition is far enoughso that its influence on the solution is negligible. The solution of Da = 140 , which has thewildest “plateau” according to Eq. (64), is analyzed at t = 40 with h = 0 . and various L . The variation of the derivative of φ at the left boundary with respect to L is shown in Fig. 3(a). When L = 15 , dφ/dx ( x = − L ) = 5 . × − . As a result, L = 15 is considered as theconverged domain length. We further investigate the convergence with respect to the grid size h . While keeping t = 40 and L = 15 , Da = 20 is selected for the mesh convergence study.Because the solution of Da = 20 has the smallest width of the diffuse interface according to Eq. (64), which leads to the highest gradient of φ among all the cases and needs the finest meshto resolve the gradient effects. By considering the solution at h = 0 . as the reference,the error e is plotted in Fig. 3 (b). It shows that our implementation is spatially second orderaccurate. The relative L error reduces to . × − when h = 0 . , which is consideredto be converged. To summarize, t = 40 , L = 15 and h = 0 . are taken as the converged parameters for the numerical simulation used in the verification.After establishing the convergence of our phase-field solver, we assess the accuracy againstEq. (64) and numerical results from [49]. The plateau-like solution at Da = 100 and thechange of h , v , w with respect to Da are shown in Fig. 4 (a)-(d), respectively. The comparisonsclearly show excellent agreements of our numerical results against previously reported ana- lytical and numerical data. The accuracy of our fully implicit variationally discretized solverfor the bistable steady convection-reaction-diffusion system in a stretching flow is successfullydemonstrated. Following the verification of our Allen-Cahn phase-field solver, we now turn our attention to the solution of Eq. (62), which can be considered as another scenario of parameter speci-fications of Eq. (63) describing the convective distortion of the diffuse interface. To form a20 a) (b)Figure 3: Convergence study of phase-field solver for a generic bistable convection-diffusion-reaction system: (a)the variation of the derivative of φ with respect to domain length L , and (b) the relative L error ( e ) as a functionof grid size h . constant ξ used in Eq. (62), the velocity field and the mobility coefficient are explicitly pre-scribed without solving the Navier-Stokes equations and the time-dependent mobility model.In the numerical simulation, the one-dimensional computational domain is taken as ˜ x ∈ [ − L, L ] , where ˜ x = x/ε is the non-dimensional coordinate. The interface thickness parameteris set to be ε = 1 . The computational domain is discretized by uniform mesh of grid size h = ∆˜ x . A zero flux boundary condition is imposed for the order parameter φ on the left andthe right boundaries. The planar interface is initialized as φ (˜ x ) = tanh(˜ x/ √ . The mobilitycoefficient of the Allen-Cahn equation is chosen as γ = 1 . The velocity is prescribed as u (˜ x ) = a ˜ x , where a ∈ [ − . , . is a constant selected according to the desired convective distortionparameter. From the problem setup, the convective distortion parameter can be calculated as ξ = a . The time step is taken as ∆ t = 0 . . The problem setup with the illustration of anextensional velocity field is shown in Fig. 5.We carry out a convergence study to minimize the discretization error so that the effect of ξ can be accurately demonstrated. To ensure that the discretization error is negligible, theconvergence of the final time when the steady-state is achieved, the length of the computationaldomain and the mesh resolution at the diffuse interface are studied with the tolerance of − for the linear GMRES and the nonlinear Newton solvers. The final time is examined for thevalues of ξ = − . and ξ = 0 . , at which the steady state solution deviates the most from the initial condition among all the cases. The grid size h = 0 . and the size of the computationaldomain L = 10 are used. By considering t = 20 as the reference, we check the relative L normof the solution at t = 18 . For ξ = − . and ξ = 0 . , the errors reduce to e = 3 . × − and e = 7 . × − , respectively. Hence the solution at t = 20 is considered to be fully-converged.We proceed to the investigation of L to check whether the zero flux boundary condition is far enough so that its influence on the solution is negligible. The solution at ξ = 0 . , which leadsto the maximum extensional distortion and the widest diffuse interface among all the cases, isanalyzed at t = 20 with h = 0 . . The derivative of the order parameter on the right boundary dφ/d ˜ x (˜ x = L ) is plotted as a function of L in Fig. 6 (a). When L = 12 , the derivative reducesto dφ/dx (˜ x = L ) = 8 . × − , which indicates an error of φ at the order of − (with h = 0 . ). Therefore it is considered as the converged domain length. The convergence of21 a) (b)(c) (d)Figure 4: Accuracy assessment of fully-implicit finite element formulation for a generic bistable convection-diffusion-reaction system: (a) steady state plateau-like solution as a function of distance x , and (b) the heightof the solution f , (c) half-width of the solution v , (d) inverse width of the diffuse interface w as a function ofDamk¨ohler number Da . the mesh resolution in the diffuse interface region is studied at t = 20 with L = 12 . Thesolution at ξ = − . is investigated, which results in the maximum compressional distortionand the highest gradient among all the cases requiring the finest mesh to resolve. By takingthe solution at ε/h = 1600 as the reference, the relative L error as a function of the mesh resolution ε/h is plotted in Fig. 6 (b). The plot confirms the second-order spatial accuracy ofour implementation. When ε/h = 800 , the error reduces to e = 6 . × − . The resolutionis deemed to be converged. To summarize, t = 20 , L = 12 and ε/h = 800 are employed tominimize the numerical error in the investigation of the effect of ξ .With the converged numerical parameters, we first consider the effect of the convectivedistortion parameter on the interface profile. As shown in Fig. 7, when ξ = 0 , the interfaceprofile from the numerical simulation tends to the hyperbolic tangent profile. When ξ > , anextensional distortion is observed. The extensional distortion increases with an increase in ξ .Compressional distortion is noted when ξ < whereby the compressional distortion increaseswith the decrease in ξ . To quantify the deviation of the interface profile from the hyperbolic22 ( φ = −
1) Ω ( φ = 1) Diffuseinterface u (˜ x ) ∼ O ( ε )2 L∂φ∂ ˜ x = 0 ∂φ∂ ˜ x = 0 Figure 5: Schematic diagram showing the computational domain for the convection of a one-dimensional planerinterface in a prescribed extensional velocity field. Ω and Ω are domains of the two phases. A zero flux boundarycondition for the order parameter is applied on the left and the right boundaries. tangent profile, we define the relative interface thickness error as: e ε = (cid:12)(cid:12)(cid:12)(cid:12) ˜ ε d − ˜ ε eq ˜ ε eq (cid:12)(cid:12)(cid:12)(cid:12) , (66)where ˜ ε d and ˜ ε eq denote the non-dimensionalized distance with respect to ε from φ = − . to φ = 0 . of the distorted interface, and of the hyperbolic tangent profile respectively. Thecoordinates of φ = 0 . and φ = − . are linearly interpolated from the numerical solution.Therelative interface thickness error as a function of ξ is shown in Fig. 10 (a).Furthermore, we examine the surface tension force calculation error due to the convectivedistortion. In the CSF model, the singular surface tension force at the interface is distributed tothe diffuse interface region by a Dirac delta function at the interface. In the current model, thefunction is given by δ S = α sf ε |∇ φ | , which should satisfy: Z ∞−∞ α sf ε (cid:18) ∂φ∂n (cid:19) dn = 1 , (67)where α sf is a constant parameter. For the hyperbolic tangent profile φ = tanh( n/ √ ε ) , theconstant parameter can be calculated as α sf = 3 √ / . When the interface profile is affected bythe convective distortion, its deviation from the hyperbolic tangent profile leads to the violationof Eq. (67). This gives rise to a relative surface tension force error quantified as: e σ = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)Z ∞−∞ α sf (cid:18) ∂φ∂ ˜ n (cid:19) d ˜ n − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (68)The error is evaluated numerically in the current study. The derivative is calculated by the cen-tral difference method at the interior points and the forward or backward difference technique at the boundary points of the computational domain. The integral is calculated by the standardtrapezoidal rule. The surface tension force error as a function of ξ is plotted in Fig. 10 (b).23 a) (b)Figure 6: Convergence of a convecting planar diffuse interface with a prescribed velocity: (a) the variation of thederivative of φ with respect to domain length L , and (b) the relative L error ( e ) as a function of mesh resolution ε/h . To demonstrate the extensional and compressional distortion in an incompressible fluidflow, we consider the convection of a curved interface in a prescribed divergence-free 2D ve-locity field. The computational domain is in the shape of an arch with r × θ ∈ [ R , R ] × [0 , π/ .A zero flux Neumann boundary condition is imposed for φ on all the boundaries. The curvedinterface is initialized as a circular arc centered at r = 0 with radius R : φ ( x, y,
0) = − tanh R − p ( x + y ) √ ε ! . (69)The curved interface is convected by a divergence-free velocity field u ( x, y ) . The horizontaland vertical components of the velocity are given by u ( x, y ) = x and v ( x, y ) = − y respectively. The mobility coefficient in the range of γ = b ∈ [4 , are investigated, where b is a constantparameter. From the problem setup, the convective distortion parameter can be calculated as ξ ( r,
0) = 1 /b at the bottom boundary and ξ ( r, π/
2) = − /b at the left boundary. In the presentcases, we use R = 0 . , R = 0 . and R = 0 . . A structured mesh of grid size ∆ r =6 . × − and ∆ θ = π/ is used for the spatial discretization. The interface thickness parameter is chosen as ε = 6 . × − . The mesh resolution can be calculated as ε/h = 10 ,which leads to a relative L convergence error at the order of − . The computational domainand the initial condition are illustrated in Fig. 8. The time step size is taken as ∆ t = 1 . × − .The numerical solutions at t = 0 . are analyzed, when the interface thickness at θ = 0 and θ = π/ reaches a constant value. The volume-preserving mean curvature flow is another source of error in the convection of acurved interface. To isolate the effect of the convective distortion from the effect of the volume-preserving mean curvature flow, the latter is minimized utilizing a small interface thicknessparameter ε . From Eq. (4), the maximum v ( x , t ) occurs at γ = 100 , t = 0 . , which can becalculated as max x ∈ Γ φI | v ( x , t ) | = 0 . . It is negligible compared to the convective velocity in the computational domain which is O (0 . . On the other hand, the minimum radius of curvature at24 igure 7: Convection of a planer interface in a prescribed velocity field: the interface profiles corresponding tovarious distortion parameter ξ . t = 0 . can be calculated as min x ∈ Γ φI R = 0 . . Since ε ≪ R , the normal velocity gradient in thenormal direction introduced by the volume-preserving mean curvature flow being O (( ε/R ) ) (by taking the derivative of v ( x , t ) with respect to R ) is negligible. Consequently, the effect ofthe volume-preserving mean curvature flow is negligible in these cases. The case at ξ = 0 . is used to illustrate the simulation results. The contour of the orderparameter from t = 0 to t = 0 . are shown in Figs. 9 (a)-(d). The time history of ˜ ε d on θ = 0 and θ = π/ , which are the bottom boundary and the left boundary respectively, are shown inFig. 9 (e). It can be observed that due to the convective distortion, the interface is expandedat θ = 0 and compressed at θ = π/ . The distortion increases during the convection of the curved interface. On the other hand, as the interface deviates further from the equilibriuminterface profile, the effect of the free energy minimization increases. When the free energyminimization opposes the convective distortion, the deviation will stop increasing, which leadsto a constant ˜ ε d . Following similar definitions of the errors in Eqs. (66), (68) and similarcalculation techniques, the relative interface thickness error and the relative surface tension force error are calculated at θ = 0 and θ = π/ . The results are shown in Fig. 10 (a)-(b)respectively. In Fig. (10), the relative interface thickness error and the relative surface tension force errorare plotted as a function of ξ for the convection of the planar and the curved interfaces. It isworth mentioning that the errors vary continuously with respect to ξ , while the singularitiesat ξ = 0 is merely a result of taking the absolute value of the errors. When the effect of thevolume-preserving mean curvature flow is minimized by reducing ε/R in the curved cases, thedependency of the errors on the convective distortion parameter ξ is almost identical for boththe cases. It shows that the relationship is consistent in one and two dimensions with differentproblem setup. Considering the numerical results and the analysis of the interface-preserving25 R R Ω ( φ = 1)Ω ( φ = − u ( x, y ) ∂φ∂n = 0 Figure 8: Schematic diagram showing the computational domain for the convection of a curved interface in aprescribed incompressible velocity field illustrated with streamlines. Ω and Ω are domains of the two phases. Azero flux Neumann boundary conditions for the order parameter is applied on all the boundaries. condition (29), we can further infer that the relationship between the errors and a constant ξ isgeneral. As shown in the figure, the errors decrease with the decrease in | ξ | . For the currentcases with a constant ξ , at the discussed location, the RMS convective distortion parameter η = | ξ | . Hence a small η can be used to reduce the error due to the convective distortion andto improve the interface-preserving capability. In Fig. 10, the interface errors e ε and e σ can befitted as functions which are proportional to η with the corresponding coefficients k ε and k σ : e ε = k ε η, e σ = k σ η, (70)where the coefficients k ε = 0 . and k σ = 0 . . The simple correlations are shownin Fig. 11. When diffuse interface is subjected to complex motions, there can be the loss of accuracy or the unresolved gradients due to the extensional and compressional distortions. Inthat case, there is a need to minimize the interface thickness error e ε . Furthermore, the controlof e σ is critical when the surface tension force and the interface dynamics play an importantrole. The RMS convective distortion parameter η can be selected as: η = min n e ε k ε , e σ k σ o for adesired value of the relative interface thickness and surface tension force errors.
5. Rising bubble problem
After analyzing the selection of the RMS convective distortion parameter, we assess theeffectiveness of our Navier-Stokes Allen-Cahn system with the time-dependent mobility modelby simulating 2D and 3D rising bubble benchmark cases. A systematic convergence study isperformed for the volume-preserving mean curvature flow and the interface-preserving capa- bility. We compare the accuracy of the diffuse interface formulation with the correspondingsharp interface method. We will show that the interface-preserving capability is necessary toguarantee an accurate solution from the diffuse interface formulation.
In this subsection, we present a well-known two-dimensional rising bubble benchmark casefor assessing the role of the time-dependent mobility model on the accurate surface tensiondynamics. The benchmark case considers the rising and deforming of an initially circular26 a) (b)(c) (d)(e)Figure 9: Convection of a curved interface with ξ = 0 . : the contours of φ at (a) t = 0 , (b) t = 0 . , (c) t = 0 . ,(d) t = 0 . , and (e) the time history of the non-dimensional interface thickness ˜ ε d on the bottom boundary θ = 0 and the left boundary θ = π/ with the comparison to the equilibrium interface thickness ˜ ε eq . a) (b)Figure 10: Dependence of the interface errors on the convective distortion parameter ξ : (a) the relative interfacethickness error e ε , and (b) the relative surface tension force error e σ . bubble immersed in a quiescent fluid. A rectangular computational domain [0 , × [0 , isconsidered for the benchmark case. The initial condition for the circular bubble is given by: φ ( x, y,
0) = − tanh R − p ( x − x c ) + ( y − y c ) √ ε ! , (71)where R = 0 . is the radius of the bubble with its center at ( x c , y c ) = (0 . , . . The quiescent fluid is initialized with zero velocity and pressure. The slip boundary condition is satisfied onthe left and the right boundaries, while the no-slip boundary condition is prescribed on the topand the bottom boundaries. A zero flux Neumann boundary condition is imposed on all theboundaries for the order parameter. The density and the dynamic viscosity of the fluid andthe bubble are selected as ρ = 1000 , ρ = 100 and µ = 10 , µ = 1 . The surface tension coefficient is set to be σ = 24 . . The gravitational acceleration is taken as g = (0 , − . .The problem setup is illustrated in Fig. 12 (a). The benchmark problem has been studied byseveral research groups employing various numerical techniques in [50]. In the current study,we consider the data from the first group in [50] for comparison purposes, which employs thefinite element method for spatial discretization and the sharp level-set method for the interface capturing. Because the problem has a symmetric axis x = 0 . , we conduct the simulation withthe right half of the computational domain. A symmetric boundary condition is imposed on theaxis of symmetry. The contour of the order parameter in the computational domain at t = 0 isshown in Fig. 12 (b). The computational domain is discretized with a uniform structured meshof grid size ∆ x = ∆ y = h . The grid size is selected according to the mesh resolution at the interface as ε/h = 1 . The time step size is chosen as ∆ t = 0 . .To quantify the mass conservation and the rising bubble dynamics, the mass of the orderparameter m , the circularity of the bubble /c , the rise velocity of the bubble V b and the center of28 a) (b)Figure 11: Correlations of the interface errors as a function of RMS convective distortion parameter η : (a) therelative interface thickness error e ε , and (b) the relative surface tension force error e σ . mass of the bubble Y b are defined as follows: m = Z Ω φ dΩ ,/c = P a /P b ,V b = R Ω v dΩ R Ω ,Y b = R Ω y dΩ R Ω , where P a is the perimeter of the circle which has the same area as the deformed bubble, P b denotes the perimeter of the bubble, v is the velocity in the Y direction and y is the Y coordinate.The defined variables and the bubble shape at t = 3 are compared in the convergence study. The volume-preserving mean curvature flow velocity disturbs the convection of the two-phase interface according to the fluid flow velocity. Therefore we perform a convergence studyto make sure that the disturbance is negligible. With the time-dependent mobility, the velocityis given by: v ( x , t ) = γ ( t ) ε κ ( x , t ) − | Γ φI ( t ) | Z Γ φI κ ( x , t ) ds ! n φL ( x , t )= 1 η F ( | ζ ( x , t ) | ) ε κ ( x , t ) − | Γ φI ( t ) | Z Γ φI κ ( x , t ) ds ! n φL ( x , t ) , x ∈ Γ φI ( t ) . (72)We refer to the quantity γ ( t ) ε as the scaling factor of the volume-preserving mean curvatureflow. According to Eq. (72), the scaling factor is affected by the user-defined parameters η ( ρ , µ ) Ω ( ρ , µ ) R Y X 1 2 u = v = 0 u = 0 u = 0 u = v = 0 (a) X Y (b)Figure 12: Two-dimensional rising bubble problem: (a) schematic diagram of the computational domain, and (b)contour of the order parameter at t = 0 . and ε . In the present convergence study, we decrease ε by a factor of from ε = 0 . to ε = 0 . , while keeping η = 0 . . Thus, the volume-preserving mean curvature flow isreduced while the RMS convective distortion parameter is kept the same. The resulting scaling factor γ ( t ) ε at t = 3 at various ε is summarized in Table 1. The circularity /c , the bubble shape φ = 0 at t = 3 , the rise velocity V b , the center of mass Y c , and their comparison with the datain [50] are shown in Fig. 13.To further quantify the results, we define the mass conservation error and the convergenceerrors as follows: e m = (cid:12)(cid:12)(cid:12)(cid:12) m t =0 − m t =3 m t =0 (cid:12)(cid:12)(cid:12)(cid:12) ,e /c = || /c − /c ref || || /c ref || ,e V b = || V b − V b, ref || || V b, ref || ,e Y b = || Y b − Y b, ref || || Y b, ref || , where |·| denotes the absolute value and ||·|| denotes the L norm of the vector. The simulationresults from the case with ε = 0 . are taken as the reference. The errors including the rel- ative interface thickness and surface tension force errors estimated by Eq. (70) are summarizedin Table 1.It can be observed from Table 1 that the mass conservation error is below and it de-creases with the reduction of ε . The total mass is well conserved. The decrease in the scalingfactor γ ( t ) ε at t = 3 along with the reduction in ε reflects the diminishing of the volume- C i r c u l a r it y (a) (b) R i s e v e l o c it y (c) C e n t e r o f m a ss (d)Figure 13: Volume-preserving mean curvature flow convergence study for 2D rising bubble benchmark case: (a)circularity of the bubble, (b) interface shape at t = 3 , (c) rise velocity, and (d) center of mass. preserving mean curvature flow. Since the flow directly affects the topology of the bubble, itsinfluence is shown clearly in the convergence of the circularity and the bubble shape. As shownin Fig. 13 (a), the circularity is reduced with the reduction of ε . This indicates that the bubbleexperiences a relatively larger deformation and a higher curvature due to the reduction of thevolume-preserving mean curvature flow. This is consistently observed in the convergence of the bubble shape shown in Fig. 13 (b). However, the bubble shape deviates from the data in[50]. The deviation can be observed in the rise velocity plotted in Fig. 13 (c) as well. This de-viation is due to the error in the surface tension force calculation introduced by the convectivedistortion. To conclude, the decreasing of the volume-preserving mean curvature flow with thedecreasing of ε cannot guarantee the accuracy of the solution, which leads to the investigation of the interface-preserving property. Following the convergence study of the volume-preserving mean curvature flow, we studythe convergence of the interface-preserving capability. From the discussion in Section 4, weknow that by decreasing the RMS convective distortion parameter η , the convective distor- able 1: Quantification of the errors for 2D rising bubble case: convergence with respect to ε at a constant η . η ε e ε e σ γ ( t ) ε ( t = 3) e m e /c e V b e Y b . .
01 0 .
074 0 .
064 1 . × − . .
005 0 .
074 0 .
064 3 . × − . . .
074 0 .
064 9 . × − . . .
074 0 .
064 2 . × − η will increase thevolume-preserving mean curvature flow. Thus we need to decrease both η and ε to improvethe interface-preserving capability and decrease the volume-preserving mean curvature flowsimultaneously. Considering this, we conduct studies with the following parameter combina-tions: η = 0 . , ε = 0 . ; η = 0 . , ε = 0 . and η = 0 . , ε = 0 . . The scaling factor of the volume-preserving mean curvature flow γ ( t ) ε at t = 3 for each combination isshown in Table 2. The circularity /c , the bubble shape at t = 3 , the rise velocity V b , the centerof mass of the bubble Y b and the comparison with the data in [50] are shown in Fig. 14. Tak-ing η = 0 . , ε = 0 . as the reference case, the errors including the relative interfacethickness and surface tension force errors are shown in Table 2. It can be observed from Table 2 that the mass is well conserved. The scaling factor ofthe volume-preserving mean curvature flow decreases with the decreasing of ε and η in the testcases. The decreasing of the RMS convective distortion parameter η reduces the errors from theconvective distortion and improves the interface preservation property. As shown in Fig. 14,with the decrease in the volume-preserving mean curvature flow and the enforcement of the interface preservation capability, the simulation results converge and match well with the datain [50]. Table 2: Quantification of the errors for 2D rising bubble case: variation of ε and η η ε e ε e σ γ ( t ) ε ( t = 3) e m e /c e V b e Y b . .
005 0 .
074 0 .
064 3 . × − . .
05 0 . .
037 0 .
032 1 . × − .
025 0 . .
018 0 .
016 1 . × − . - - - We further emphasize the importance of the interface-preserving capability by comparingthe results of the 2D rising bubble benchmark case simulated by the constant mobility model with γ = 1 and the time-dependent mobility with η = 0 . . The complete computational domain [0 , × [0 , is used in the simulations. The interface thickness parameter is set to be ε = 0 . and a uniform structured mesh ∆ x = ∆ y = h of grid size ε/h = 1 is employed for the spatialdiscretization via linear finite elements. The time step size is taken as ∆ t = 0 . .The contour of the order parameter at t = 3 superimposed on the mesh simulated by the constant and the time-dependent mobility are shown in Figs. 15 (a) and (b), respectively. Asobserved in Fig. 15 (a), when γ = 1 , the interface-preserving capability is not sufficient to keepthe interface profile against the convective distortion. Therefore, at the bottom of the bubble,the interface is subjected to an observable extensional distortion, which leads to an excessivelylow Laplace pressure. This changes the shape of the bubble, decreases the buoyancy force, i r c u l a r it y (a) (b) R i s e v e l o c it y (c) C e n t e r o f m a ss (d)Figure 14: The interface-preserving capability convergence study for 2D rising bubble benchmark case: (a) circu-larity of the bubble (b) interface shape at t = 3 (c) rise velocity, and (d) center of mass. and further reduces the rise velocity as shown in Fig. 15 (c). On the contrary, when the time-dependent mobility model with η = 0 . is used, the interface profile is preserved well asshown in Fig. 15 (b), which gives a correct bubble shape and the surface tension force. Tofurther justify the above statements, we quantify the convective distortion by calculating theRMS convective distortion parameter as η = F ( | ζ ( x , t ) | ) /γ . As shown in Fig .15 (d), a larger η representing insufficient interface-preserving capability is observed in the simulation withthe constant mobility compared to the time-dependent mobility. The comparison shows thatthe proposed time-dependent model provides an approach to estimate as well as control theconvective distortion. The 3D rising bubble benchmark case is a generalization of the 2D rising bubble case withincreasing complexity and practicality due to the 3D topology and motion of the bubble. Weuse the benchmark case to further assess the Navier-Stokes Allen-Cahn CSF system with thetime-dependent mobility model. The benchmark case considers the rising and deforming of an33 a) (b) R i s e v e l o c it y (c) Constant mobility modelTime-dependent mobility model (d)Figure 15: Comparison of the constant and the time-dependent mobility model in 2D rising bubble case: contourof the order parameter at t = 3 simulated with (a) constant mobility coefficient γ = 1 , (b) time-dependent mobilitymodel at η = 0 . , and the difference in (c) rise velocity, and (d) the time history of the RMS convective distortionparameter η . Ω ( ρ , µ ) Ω ( ρ , µ ) R (a) (b)Figure 16: Three-dimensional rising bubble problem: (a) schematic diagram of the computational domain, and(b) contour of the order parameter at t = 0 . initially spherical bubble in a cuboid tank occupying the spatial domain Ω ∈ [0 , × [0 , × [0 , . The phase-field function describing the bubble is initialized as: φ ( x, y, z,
0) = − tanh R − p ( x − x c ) + ( y − y c ) + ( z − z c ) √ ε ! , (73)where R = 0 . is the radius of the bubble with its center at ( x c , y c , z c ) = (0 . , . , . . Theno-slip boundary condition and the zero flux Neumann boundary condition are imposed onall the boundaries for the velocity and the order parameter respectively. The density and theviscosity of the fluid and the bubble are taken as ρ = 1000 , ρ = 100 , µ = 10 , µ = 1 . Thesurface tension coefficient is taken as σ = 24 . . The gravitational acceleration is set to be g = (0 , − . , . The problem setup is illustrated in Fig. 16 (a). The 3D rising bubble benchmarkcase is investigated in [51] by several research groups. We consider the data from the secondgroup in [51] for comparison purposes, in which the finite difference method and the sharplevel-set method are employed for the discretization and the interface capturing respectively.As the case is symmetric with respect to planes x = 0 . and z = 0 . , we simulate one quarter of the computational domain, as shown in Fig. 16 (b) with the symmetric boundary conditionimposed on the symmetric planes. The computational domain is discretized with a uniformstructured mesh of grid size ∆ x = ∆ y = ∆ z = h . The mesh resolution at the interface isselected as ε/h = 1 , while the time step is taken as ∆ t = 0 . .We define the following variables to assess the simulation results quantitatively: the totalmass of the order parameter m, the sphericity of the bubble /s , the diameter of the bubble in Xdirection D x and Y direction D y , the rise velocity of the bubble V b and the center of mass of35he bubble Y b , which are given by: m = Z Ω φd Ω ,/s = A a /A b ,D x = max( { x | x ∈ Ω } ) − min( { x | x ∈ Ω } ) ,D y = max( { y | y ∈ Ω } ) − min( { y | y ∈ Ω } ) ,V b = R Ω vd Ω R Ω d Ω ,Y b = R Ω yd Ω R Ω d Ω , where A a is the area of the sphere which has the same volume as the deformed bubble, A b denotes the surface area of the bubble, x and y are the coordinates in X and Y directionsrespectively, and v is the velocity in the Y direction. The defined variables are compared in theconvergence study.Similar to the convergence study of the 2D rising bubble case, we decrease η and ε to im-prove the interface preservation capability and decrease the volume-preserving mean curvature flow simultaneously.The following combinations are tested: η = 0 . , ε = 0 . ; η = 0 . , ε =0 . and η = 0 . , ε = 0 . . The scaling factor of the volume-preserving mean curvatureflow γ ( t ) ε at t = 3 for each combination is tabulated in Table 3. The defined variables andtheir comparison to the data in [51] are shown in Fig. 17.To further quantify the results, we define the mass conservation error and convergence errorsas: e m = (cid:12)(cid:12)(cid:12)(cid:12) m t =0 − m t =3 m t =0 (cid:12)(cid:12)(cid:12)(cid:12) ,e /s = || /s − /s ref || || /s ref || ,e D x = || D x − D x, ref || || D x, ref || ,e D y = || D y − D y, ref || || D y, ref || ,e V b = || V b − V b, ref || || V b, ref || ,e Y b = || Y b − Y b, ref || || Y b, ref || . The simulation results of η = 0 . , ε = 0 . are taken as the reference. The errors including the relative interface thickness and surface tension force errors are tabulated in Table 3.It can be observed that the total mass is well conserved with less than relative error.The relative interface thickness and surface tension force errors decreases with the decreasein the RMS convective distortion parameter η . The reduction of the volume-preserving meancurvature flow is reflected by the reduction of its scaling factor γ ( t ) ε at t = 3 in current cases. As a result, lower sphericity representing larger deformation and higher curvature is36 able 3: Quantification of the errors for the 3D rising bubble benchmark η ε e ε e σ γ ( t ) ε ( t = 3) e m e /s e D x e D y e V b e Y b . .
01 0 .
074 0 .
064 1 . × − . .
05 0 .
005 0 .
037 0 .
032 8 . × − . .
025 0 . .
018 0 .
016 5 . × − . - - - - -observed in Fig. 17 (a). This is further confirmed in Fig. 17 (b), where a larger D x and asmaller D y representing more deviation from the spherical shape are observed. With both theimprovement of the interface-preserving capability and the decrease in the volume-preservingmean curvature flow, the converged simulation results agree well with the data in [51]. S ph e r i c it y (a) D i a m e t e r (b) R i s e v e l o c it y (c) C e n t e r o f m a ss (d)Figure 17: Convergence study for the 3D bubble rising problem: (a) sphericity of the bubble, (b) bubble diametersin X and Y directions, (c) rise velocity, and (d) center of mass.
6. Two rising bubbles merging with a free surface
In this section, we demonstrate the applicability of the proposed model in the case of tworising bubbles merging with a free surface, in which complicated topological changes of the37nterface and dynamics of the bubble-bubble and bubble-free surface interaction occur in anunstructured finite element mesh. The case considers the rising of two vertically aligned spher-ical bubbles driven by gravitational force and the merging of bubbles with the free surfacein a cuboid tank. The computational domain is taken as Ω ∈ [0 , × [0 , × [0 , . Theno-slip boundary condition and the zero flux Neumann boundary condition are applied on allthe boundaries for the velocity and the order parameter, respectively. The order parameter isinitialized for the bubbles and the free surface as: φ ( x, y, z,
0) = − tanh R l − q ( x − x lc ) + ( y − y lc ) + ( z − z lc ) √ ε − tanh R u − q ( x − x uc ) + ( y − y uc ) + ( z − z uc ) √ ε − tanh (cid:18) y − y wl √ ε (cid:19) − , (74)where R l = 0 . is the radius of the lower bubble with its center at ( x lc , y lc , z lc ) = (0 . , . , . , R u = 0 . is the radius of the upper bubble with its center at ( x uc , y uc , z uc ) = (0 . , , . , y wl = 1 . is the water level of the free surface. The density and viscosity of the fluid and thebubbles are taken as ρ = 1000 , ρ = 100 , µ = 10 , µ = 1 . The surface tension coefficient is chosen as σ = 24 . . The gravitational acceleration is set to be g = (0 , − . , . The problemdefinition is illustrated in Fig. 18 (a). The interface thickness parameter is selected as ε = 0 . to reduce the volume-preserving mean curvature flow. In the time-dependent mobility model,the RMS convective distortion parameter η = 0 . is used to get an accurate surface tensionforce calculation. The above combinations of ε and γ has been proven to be accurate in the sim- ulation of the 3D rising bubble case, which uses the same physical parameters. A non-uniformunstructured mesh consisting of 7,077,043 nodes and 45,078,392 tetrahedrons is employed forthe spatial discretization. The mesh at the plane x = 0 . is shown in Fig. 18 (b). The time stepsize is selected as ∆ t = 0 . . The evolution of the interface φ = 0 is shown in Fig. 19. Thecomplex topological changes including the rising of the bubbles, the merging of the bubbles with the free surface and the wave formation at the free surface in the merging process can beobserved clearly using our 3D phase-field Navier-Stokes solver with the unstructured mesh.
7. Conclusion
In the present work, a variational interface-preserving Allen-Cahn phase-field formulationrelying on a novel time-dependent mobility model has been developed for accurate surface ten- sion force calculation. By writing the convective Allen-Cahn equation in a non-dimensionalmoving orthogonal curvilinear coordinate system, we have derived the governing equation forthe interface profile. We have identified the convective distortion term and the effective parame-ter determining the deviation of the diffuse interface profile from the hyperbolic tangent profilein the governing equation. A time-dependent mobility model to control the convective distor- tion parameter in the diffuse interface region and to preserve the hyperbolic tangent profile hasbeen constructed accordingly. Following the verification of our implicit PPV-based steady-stateAllen-Cahn solver for a generic bistable convection-reaction-diffusion system, we established38 ( ρ , µ ) Ω ( ρ , µ ) R l R u XZY Ω Ω y wl (a) (b)Figure 18: Two rising bubbles merging with a free surface problem: (a) schematic diagram showing the computa-tional domain, and (b) unstructured finite element mesh at the plane x = 0 . . a) (b) (c)(d) (e) (f)(g) (h) (i)Figure 19: Two rising bubbles merging with a free surface: the evolution of the interface φ = 0 at different timeinstants t = (a) 0.1, (b) 1.1, (c) 1.6, (d) 1.7, (e) 1.8, (f)1.9, (g) 2.2, (h) 2.4, and (i) 2.7. the solutions of the proposed model in two- and three-dimensional rising bubble benchmarkcases against the sharp interface counterparts. Through the assessment, it has been shownthat the interface preservation achieved by the proposed model and the minimization of thevolume-preserving mean curvature flow realized by decreasing the interface thickness param-eter are essential for the accurate surface tension dynamics. Finally, by simulating two rising bubbles merging with a free surface, we have shown that the proposed technique is applicablein a practical problem involving complex topology changes of the interface in an unstructuredmesh and complex dynamics involving bubble-bubble and bubble-free surface interactions. Acknowledgement
The authors would like to acknowledge the Natural Sciences and Engineering Research
Council of Canada (NSERC) for the funding. This research was supported in part throughcomputational resources and services provided by Advanced Research Computing at the Uni-versity of British Columbia.
Appendix A. Quantification of the convective distortion intensity
In this appendix, we demonstrate the quantification of the intensity of the convective dis- tortion by the normal velocity gradient in the normal direction ζ . Consider two level sets of φ = φ and φ = φ . Due to the finite thickness of the region, the convective velocity of thetwo level sets can be different. The tangential velocities are ignored since they do not affectthe propagation of level sets. The normal velocity of the level sets are denoted as u n ( x , t ) and u n ( x , t ) , where x and x are in the same normal axis with identical tangential coordinates. Suppose φ is continuous and differentiable with respect to x ∈ ( x , x ) , we want to quantifythe convective distortion intensity at x .Consider the initial distance between the two level sets at the discussed normal axis as ε i = || x − x || . After a small time increment ∆ t , the thickness is distorted due to thedifference in the normal velocity, which can be linearly approximated as: ε e = ε i + ( u n ( x , t ) − u n ( x , t ))∆ t, (A.1)where ε e represents the distance between the two level sets at the end of the time increment, asshown in Fig. A.1. The intensity of the convective distortion at x ∈ ( x , x ) can be consideredas the relative change of the thickness per unit time due to the difference in the convectivevelocity: lim ε i → t → ε e − ε i ε i ∆ t = lim ε i → t → u n ( x , t ) − u n ( x , t ) ε i = ∂u n ∂n ( x , t ) = ζ ( x , t ) (A.2) Appendix B. Frame independent form of time-dependent mobility
The normal velocity gradient in the normal direction ζ ( x , t ) in Eq. (31) can be written in aframe independent form: ζ ( x , t ) = ∇ ( u · n φL ) · n φL . (B.1)41 n ( x , t ) u n ( x , t ) φ = φ φ = φ ε i ∆ t u n ( x , t ) u n ( x , t ) φ = φ φ = φ ε e Figure A.1: Schematic diagram of the convective distortion of the thickness between the level sets of φ = φ and φ = φ due to different normal velocities. To simplify the equation, we expand the normal velocity gradient in the normal direction as: ∇ ( u · n φL ) · n φL = n φL · ∇ u · n φL + u · ∇ n φL · n φL + u × ( ∇ × n φL ) · n φL + n φL × ( ∇ × u ) · n φL where the second, the third and the fourth terms are zero. For example, the second term: u · ∇ n φL · n φL = u · ( ∇ n φL · n φL ) = u · ∇ ( n φL · n φL ) = 0 and the third term ∇ × n φL = ∇ × (cid:16) ∇ φ |∇ φ | (cid:17) = ∇ (cid:16) |∇ φ | (cid:17) × ∇ φ + |∇ φ | ∇ × ∇ φ . Notice that with the assumption of uniforminterface profile along the interface, / |∇ φ | is a constant on the level sets of φ . Thus both ∇ (1 / |∇ φ | ) and ∇ φ are normal to the level sets of φ , which leads to ∇ (1 / |∇ φ | ) × ∇ φ = 0 .The curl of gradient of a scalar field is identically zero. As a result, the third term vanishes.The last term equals to zero since for arbitrary vectors a and b , ( a × b ) · a = 0 . As a result,the normal velocity gradient in the normal direction can be simplified as: ζ ( x , t ) = n φL · ∇ u · n φL . (B.2)and by substituting n φL = ∇ φ/ |∇ φ | , we get the frame-independent form of the normal velocitygradient in the normal direction: ζ ( x , t ) = ∇ φ · ∇ u · ∇ φ |∇ φ | . (B.3)Finally, the frame independent form of the time-dependent mobility model is given by: γ ( t ) = 1 η F (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ∇ φ · ∇ u · ∇ φ |∇ φ | (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) , (B.4)where F ( ϕ ( x , t )) = q R ( ϕ ( x ,t )) d Ω R d Ω , x ∈ Γ φDI ( t ) . The frame independent form facilitates itsnumerical implementation in Cartesian coordinate system. References [1] C. E. Brennen, Cavitation and bubble dynamics, Cambridge University Press, 2014.422] A. Smirnov, I. Celik, S. Shi, LES of bubble dynamics in wake flows, Computers & Fluids34 (3) (2005) 351–373.[3] B. Mallat, G. Germain, J.-Y. Billard, B. Gaurier, A 3D study of the bubble sweep-down phenomenon around a 1/30 scale ship model, European Journal of Mechanics-B/Fluids72 (2018) 471–484.[4] C. Hirt, A. Amsden, J. Cook, An arbitrary Lagrangian Eulerian computing method for allflow speeds, Journal of Computational Physics 135 (2) (1997) 203–216.[5] S. O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, Journal of Computational Physics 100 (1992) 25–37.[6] W. J. Rider, D. B. Kothe, Reconstructing volume tracking, Journal of ComputationalPhysics 141 (2) (1998) 112–152.[7] C. W. Hirt, B. D. Nichols, Volume of fluid (VOF) method for the dynamics of free bound-aries, Journal of Computational Physics 39 (1) (1981) 201–225. [8] R. Malladi, J. A. Sethian, B. C. Vemuri, Shape modeling with front propagation: A levelset approach, IEEE transactions on pattern analysis and machine intelligence 17 (2) (1995)158–175.[9] J. U. Brackbill, D. B. Kothe, C. Zemach, A continuum method for modeling surfacetension, Journal of Computational Physics 100 (2) (1992) 335–354. [10] J. D. van der Waals, The thermodynamic theory of capillarity under the hypothesis of acontinuous variation of density, Journal of Statistical Physics 20 (2) (1979) 200–244.[11] D. Jacqmin, Calculation of two-phase Navier–Stokes flows using phase-field modeling,Journal of Computational Physics 155 (1) (1999) 96–127.[12] E. Olsson, G. Kreiss, A conservative level set method for two phase flow, Journal of
Computational Physics 210 (1) (2005) 225–246.[13] O. Desjardins, V. Moureau, H. Pitsch, An accurate conservative level set/ghost fluidmethod for simulating turbulent atomization, Journal of Computational Physics 227 (18)(2008) 8395–8416.[14] D. Hartmann, M. Meinke, W. Schr¨oder, On accuracy and efficiency of constrained reini- tialization, International Journal for Numerical Methods in Fluids 63 (11) (2010) 1347–1358.[15] A. K. Tornberg, B. Enhquist, A finite element based level set method for multiphase flowapplications, Comput. Visual. Sci. 3 (2000) 93–101.[16] P. H. Chiu, Y. T. Lin, A conservative phase field method for solving incompressible two- phase flows, Journal of Computational Physics 230 (1) (2011) 185–204.[17] V. Joshi, R. K. Jaiman, A positivity preserving and conservative variational scheme forphase-field modeling of two-phase flows, Journal of Computational Physics 360 (2018)137–166. 4318] J. W. Cahn, J. E. Hilliard, Free energy of a nonuniform system. I. Interfacial free energy,
The Journal of Chemical Physics 28 (2) (1958) 258–267.[19] J. W. Cahn, On spinodal decomposition, Acta Metallurgica 9 (9) (1961) 795–801.[20] S. M. Allen, J. W. Cahn, A microscopic theory for antiphase boundary motion and itsapplication to antiphase domain coarsening, Acta Metallurgica 27 (6) (1979) 1085–1095.[21] J. Kim, S. Lee, Y. Choi, S. M. Lee, D. Jeong, Basic principles and practical applications of the Cahn–Hilliard equation, Mathematical Problems in Engineering 2016 (2016).[22] J. Rubinstein, P. Sternberg, Nonlocal reaction-diffusion equations and nucleation, IMAJournal of Applied Mathematics 48 (3) (1992) 249–264.[23] M. Brassel, E. Bretin, A modified phase field approximation for mean curvature flowwith conservation of the volume, Mathematical Methods in the Applied Sciences 10 (34) (2011) 1157–1180.[24] Y. Sun, C. Beckermann, Sharp interface tracking using the phase-field equation, Journalof Computational Physics 220 (2) (2007) 626–653.[25] Z. Chai, D. Sun, H. Wang, B. Shi, A comparative study of local and nonlocal Allen-Cahnequations with mass conservation, International Journal of Heat and Mass Transfer 122 (2018) 631–642.[26] J. Shen, X. Yang, Numerical approximations of Allen-Cahn and Cahn-Hilliard equations,Discrete Contin. Dyn. Syst 28 (4) (2010) 1669–1691.[27] L. Modica, The gradient theory of phase transitions and the minimal interface criterion,Archive for Rational Mechanics and Analysis 98 (2) (1987) 123–142. [28] J. Rubinstein, P. Sternberg, J. B. Keller, Fast reaction, slow diffusion, and curve shorten-ing, SIAM Journal on Applied Mathematics 49 (1) (1989) 116–133.[29] P. Yue, J. J. Feng, C. Liu, J. Shen, A diffuse-interface method for simulating two-phaseflows of complex fluids, Journal of Fluid Mechanics 515 (2004) 293–317.[30] Y. Li, J.-I. Choi, J. Kim, A phase-field fluid modeling and computation with interfacial profile correction term, Communications in Nonlinear Science and Numerical Simulation30 (1-3) (2016) 84–100.[31] Y. Zhang, W. Ye, A flux-corrected phase-field method for surface diffusion, Communica-tions in Computational Physics 22 (2) (2017) 422–440.[32] G. Soligo, A. Roccon, A. Soldati, Mass-conservation-improved phase field methods for turbulent multiphase flow simulation, Acta Mechanica 230 (2) (2019) 683–696.[33] V. Joshi, R. K. Jaiman, An adaptive variational procedure for the conservative and pos-itivity preserving Allen-Cahn phase-field model, Journal of Computational Physics 366(2018) 478–504. 4434] V. Joshi, R. K. Jaiman, A hybrid variational Allen-Cahn/ALE scheme for the coupled analysis of two-phase fluid-structure interaction, International Journal for NumericalMethods in Engineering 117 (4) (2019) 405–429.[35] V. Joshi, R. K. Jaiman, A positivity preserving variational method for multi-dimensionalconvection–diffusion–reaction equation, Journal of Computational Physics 339 (2017)247–284. [36] R. G. Bartle, D. R. Sherbert, Introduction to real analysis, Vol. 2, Wiley New York, 2000.[37] J. Kim, Phase-field models for multi-component fluid flows, Communications in Compu-tational Physics 12 (3) (2012) 613–661.[38] J. Kim, A continuous surface tension force formulation for diffuse-interface models, Jour-nal of Computational Physics 204 (2) (2005) 784–804. [39] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, G. Zanetti, Modelling merging andfragmentation in multiphase flows with surfer, Journal of Computational Physics 113 (1)(1994) 134–147.[40] J. Chung, G. M. Hulbert, A time integration algorithm for structural dynamics with im-proved numerical dissipation: The generalized- α method, Journal of Applied Mechanics
60 (2) (1993) 371–375.[41] T. J. Hughes, G. N. Wells, Conservation properties for the Galerkin and stabilised forms ofthe advection–diffusion and incompressible Navier–Stokes equations, Computer Methodsin Applied Mechanics and Engineering 194 (9-11) (2005) 1141–1159.[42] M. C. Hsu, Y. Bazilevs, V. M. Calo, T. E. Tezduyar, T. J. Hughes, Improving stability of stabilized and multiscale formulations in flow simulations at small time steps, ComputerMethods in Applied Mechanics and Engineering 199 (13-16) (2010) 828–840.[43] F. Shakib, T. J. Hughes, Z. Johan, A new finite element formulation for computationalfluid dynamics: X. The compressible Euler and Navier-Stokes equations, Computer Meth-ods in Applied Mechanics and Engineering 89 (1-3) (1991) 141–219. [44] A. N. Brooks, T. J. Hughes, Streamline upwind/Petrov-Galerkin formulations for con-vection dominated flows with particular emphasis on the incompressible Navier-Stokesequations, Computer Methods in Applied Mechanics and Engineering 32 (1-3) (1982)199–259.[45] I. Harari, T. J. Hughes, What are C and h?: Inequalities for the analysis and design of finite element methods, Computer Methods in Applied Mechanics and Engineering 97 (2)(1992) 157–192.[46] C. Johnson, Numerical solutions of partial differential equations by the finite elementmethod, Cambridge University Press, 1987.[47] R. Jaiman, M. Guan, T. Miyanawala, Partitioned iterative and dynamic subgrid-scale methods for freely vibrating square-section structures at subcritical reynolds number,Computers & Fluids 133 (2016) 68–89. 4548] Y. Saad, M. H. Schultz, Gmres: A generalized minimal residual algorithm for solvingnonsymmetric linear systems, SIAM Journal on Scientific and Statistical Computing 7 (3)(1986) 856–869. [49] S. M. Cox, G. A. Gottwald, A bistable reaction-diffusion system in a stretching flow,Physica D: Nonlinear Phenomena 216 (2) (2006) 307–318.[50] S. R. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, L. Tobiska,Quantitative benchmark computations of two-dimensional bubble dynamics, InternationalJournal for Numerical Methods in Fluids 60 (11) (2009) 1259–1288.750