A higher-order finite element method with unstructured anisotropic mesh adaption for two phase flows with surface tension
AA higher-order finite element method with unstructured anisotropicmesh adaption for two phase flows with surface tension
Modesar Shakoor ∗ and Chung Hae ParkIMT Lille Douai, Institut Mines-T´el´ecom, Univ. Lille, Centre for Materials and Processes,F-59000 Lille, FranceOctober 27, 2020 Abstract
A novel finite element framework is proposed for the numerical simulation of two phase flows withsurface tension. The Level-Set (LS) method with piece-wise quadratic (P2) interpolation for the liquid-gas interface is used in order to reach higher-order convergence rates in regions with smooth interface.A balanced-force implementation of the continuum surface force model is used to take into account thesurface tension and to solve static problems as accurately as possible. Given that this requires a balancebetween the discretization used for the LS function, and that used for the pressure field, an equal-orderP2/P2/P2 scheme is proposed for the Navier-Stokes and LS advection equations, which are stronglycoupled with each other. This fully implicit formulation is stabilized using the residual-based variationalmultiscale framework. In order to improve the accuracy and obtain optimal convergence rates with aminimum number of elements, an anisotropic mesh adaption method is proposed where the unstructuredmesh is kept as fine as possible close to the zero iso-value of the P2 LS function. Elements are automati-cally stretched in regions with flat interface in order to keep the complexity fixed during the simulation.The accuracy and efficiency of this approach are demonstrated for two and three dimensional simulationsof a rising bubble.
Keywords: multiphase flow, level-set, anisotropic mesh, variational multiscale, higher-order, balanced-force
Trending industrial applications such as advanced manufacturing ( e.g. , liquid composite moulding, additivemanufacturing) feature complex multiphase flows where capillary effects have a significant influence at differ-ent scales [1, 2, 3]. An increasing research effort has been made in the last decade towards the developmentof robust and reliable Finite Element (FE) simulation tools for the modeling of multiphase flows with surfacetension.First, stabilized FE methods are necessary in order to solve Navier-Stokes equations for a wide range of flowregimes, including turbulent [4]. These methods can also be used to satisfy the inf-sup condition when usingincompatible velocity-pressure discretizations [5]. Recent developments based on the Variational MultiScale(VMS) framework are particularly relevant in that regard [6]. For instance, the Residual-Based VMS (RB-VMS) framework has been proven to be applicable both to reaction-diffusion-advection equations and to theincompressible Navier-Stokes equations in convection-dominated regimes for various choices of discretization[7, 8]. ∗ Corresponding author a r X i v : . [ c s . C E ] O c t econd, as shown in Ref. [9], advanced mesh adaption and motion techniques are required to maintain aconformal FE mesh at interfaces between the different phases of the flow during the simulation. Instead,a front-tracking method is generally employed in order to represent and track these moving and deforminginterfaces. Most surface tension models additionally introduce a mean curvature dependent surface forceterm at liquid-gas interfaces [10]. The volume of fluid method [11], which tracks the local volume fractionof each phase during the simulation, is hence not sufficient, because it does not capture the morphology ofthe interface. The interface morphology can be captured by using the Level-Set (LS) method [12], whichapproximates the boundary of each phase through its signed distance function, a.k.a. LS function. Indeed,the normal vector and the mean curvature of an interface are related to the first and second derivatives of itsdistance function [13]. The properties of volume conservation in the volume of fluid method are however lostwhen switching to the LS method, which has led researchers to consider combinations of those two methods[14]. Alternatively, higher-order discretizations using discontinuous Galerkin methods [15] or isogeometricanalysis [8], anisotropic mesh adaption [16, 1] and semi-Lagrangian particles [17] have been considered todeal with the same issue, as well as to improve the order of accuracy at the interface.Because the liquid-gas interface is usually not discretized with a conformal FE mesh, the implementation ofthe surface force term is not straightforward. Both the discontinuity in material properties ( i.e. , density andviscosity) and the surface tension term can be integrated accurately by reconstructing the separate liquidand gas domains and the interface at elements intersected by the interface using the extended FE method[18, 19]. This interface reconstruction and integration step can be avoided using the Continuum SurfaceForce (CSF) model, where the surface force is transformed into a volume force distributed in a small layeraround the interface through the Dirac delta function [10]. Similarly, the discontinuity in material propertiescan be distributed in this small layer through a regularized Heaviside function [13]. It is important to pointout that the regularized Heaviside function depends on the LS function, and that the delta function dependson the first derivative of the regularized Heaviside function, which can be computed analytically [10, 13].The regularized transition of the material properties relies on the distance property of LS functions [13].This distance property is usually lost during interface motion and deformation. Thus, it should be restoredthrough an operation called LS reinitialization. This can be achieved directly through geometric distancecomputation [20] or indirectly by solving an Eikonal equation using fast-marching methods [21]. Alterna-tively, indirect LS reinitialization methods based on the solution of nonlinear hyperbolic, parabolic or ellipticequations have been proposed in the literature [22]. The choice of LS reinitialization method is not obvi-ous and there is no method that has been widely accepted in the literature. Geometric reinitialization istheoretically more robust, because the exact distance to the interface can be computed [20]. However, itis computationally expensive and its generalization to different FE types and orders is not straightforward[23]. Similarly, the generalization of fast-marching methods to unstructured FE meshes requires generalizedupwinding criteria [21]. Indirect methods based on a nonlinear hyperbolic equation rely on an iterativeprocess which has been shown to be unstable in the presence of interface singularities [20]. Recently, LSreinitialization methods based on nonlinear parabolic or elliptic equations have been developed. Even if theirrobustness has not been demonstrated in the frame of multiphase flow simulations yet [22], these indirectmethods are promising as they can easily be implemented for any FE type and order.Spurious non-physical currents near the interface have been reported in the literature, even for static cases[24, 23, 25, 26]. As analyzed in details in Refs. [25, 27], these currents originate from a mismatch betweenthe pressure discretization and the surface tension model, and numerical errors in the mean curvature calcu-lation. To address the first issue, a balanced-force algorithm can be employed. In this algorithm, the deltafunction and the normal vector in the CSF model are replaced by the numerical gradient of the regularizedHeaviside function, and this gradient as well as the pressure gradient are projected to the same discretizationspace [24]. To avoid those projections, the same discretization can be used for the LS function and the pres-sure. Since the pressure discretization space cannot be larger than that of the velocity, the balanced-forceimplementation of higher-order LS methods is more difficult.Last, interface motion in the LS method is modeled through advection equations for LS functions. Theadvection velocity is no other than the flow velocity computed by the flow solver. The flow model includesthe CSF model, which depends on the second derivatives of the LS function and, in a balanced-force im-2lementation, on the first derivatives of the LS-dependent regularized Heaviside function. Thus, a choiceof time stepping scheme and coupling should be made. To avoid severe restrictions on the time step [13],various implicit schemes have been considered in the literature. It has been shown that both a semi-implicitand a fully implicit implementation of the CSF model can be achieved without strong coupling of the flowand LS advection solvers [28]. To eliminate all restrictions on the time step, fully implicit implementationswith strong coupling of the two solvers have been proposed [29, 8].In this work, we propose an adaptive fully implicit RBVMS method for two phase flow with surface tensionrelying on a piece-wise quadratic (P2) discretization of the velocity, the pressure and the LS function as well,in order to reduce the computational cost while maintaining a good accuracy of the numerical solutions.All variables, i.e. , the velocity, the pressure and the LS function are solved simultaneously because Navier-Stokes equations with the CSF model are strongly coupled to the LS advection equation, as opposed to mostconventional approaches. The use of the same discretization space for all variables is only possible becauseof RBVMS stabilization, and enables us to obtain a balanced-force implementation in a straightforwardmanner without the procedure of pressure gradient projection. Moreover, a large time step can be adoptedthrough the fully implicit RBVMS method while avoiding divergence of the solution. In order to improvethe accuracy near the interface, the unstructured mesh is dynamically adapted during the simulation basedon an anisotropic P2 error estimator.The governing equations for incompressible two phase flow with surface tension are summarized in Sec. 2,and the proposed numerical method is presented in Sec. 3. In Sec. 4, the performance of the proposedmethod in terms of accuracy and efficiency is demonstrated for two and three dimensional simulations of arising bubble. It is shown that our fully implicit implementation leads to a reduction of the number of timesteps compared to other methods proposed in the literature, while anisotropic unstructured mesh adaptionallows us to achieve a more optimal FE discretization with a higher-order convergence rate. The computation domain is denoted as Ω ⊂ R d , d = 2 ,
3. The Navier-Stokes momentum and mass conserva-tion equations for isothermal, Newtonian, incompressible flow in Ω are ρ (cid:18) ∂ v ∂t + v . ∇ v (cid:19) −∇ .σ ( v , p ) − f g + f s = , ∇ . v = 0 . (1)In these equations v is the velocity vector, p the pressure, t the time, ρ the mass density, and µ the viscosity.The Cauchy stress tensor σ is defined as σ ( v , p ) = µ ( ∇ v + ∇ T v ) − p I , (2)where I is the identity matrix. The force f g = − ρg e y is the gravity force, with g the gravity acceleration.The force f s is the surface tension force and is described in Subsec. 2.3.Those equations will be accompanied by either Dirichlet or homogeneous Neumann boundary conditions inthe paper. The computation domain is split into two parts Ω and Ω , each occupied by a different fluid, withΩ ∪ Ω = Ω, and Ω ∩ Ω = Γ , . Both domains as well as their interface evolve in time according to the motion of the two fluids:Ω = Ω ( t ) , Ω = Ω ( t ), and Γ , = Γ , ( t ) . , , we introduce the LS function φ [12] which is the signed distancefunction to Γ , : φ ( x , t ) = + dist ( x , Γ , ( t )) , x ∈ Ω ( t ) , − dist ( x , Γ , ( t )) , x ∈ Ω ( t ) , , x ∈ Γ , ( t ) . (3)The first (resp. second) fluid hence occupies the part of the domain with positive (resp. negative) φ ,while the interface is defined as the zero isolevel of φ . The evolution of the interface can then be modeledby advecting φ with the flow velocity: ∂φ∂t + v . ∇ φ = 0 . (4)Fluid properties in Eq. (1) can now be defined in terms of the LS function [30]: ρ ( x , t ) = ρ ( φ ( x , t )) = H ( φ ( x , t )) ρ + (1 − H ( φ ( x , t ))) ρ µ ( x , t ) = µ ( φ ( x , t )) = H ( φ ( x , t )) µ + (1 − H ( φ ( x , t ))) µ (5)where ρ , ρ and µ , µ are the mass densities and viscosities of each fluid, respectively, and H is the Heavisidefunction: H ( φ ) = , φ > , , φ = 0 , , φ < . (6) The surface tension force f s can be expressed through the CSF model [10] as f s = σ s κ s n s δ (7)where σ s is the surface tension coefficient, κ s is the mean curvature of the interface, n s its normal vector,and δ = d H d φ is the so-called interface density. The computation of this derivative requires a regularizationof the Heaviside function, as described in Subsec. 3.1.The normal vector and the mean curvature can be directly computed from LS function φ using n s = n s ( φ ) = ∇ φ |∇ φ | , κ s = κ s ( φ ) = −∇ . n s ( φ ) = ∇ φ. ∇∇ φ. ∇ φ − |∇ φ | tr( ∇∇ φ ) |∇ φ | , (8)where ∇∇ φ is the Hessian matrix of φ and | . | the Euclidean norm. Due to the sign of the LS function inEq. (3), the normal vector will be pointing towards phase 1, which means that the liquid phase is phase 1while the gas phase is phase 2.The numerical formulation described in the following section solves Eqs. (1) and (4) to approximate { v , p, φ } . It can be observed that the LS advection equation depends linearly on the flow velocity, whilethe momentum equation involves a nonlinear dependence of the mass density and the viscosity on the LSfunction. This is added to the nonlinearity of the auto-advection term and of the surface tension force. In this section, we present the different numerical methods to solve the equations described in Sec. 2.Those equations are solved using the RBVMS framework where we choose equal-order FE discretizations { v h , p h , φ h } for the velocity, pressure and LS solutions. This approximation space is denoted by V h . Theparameter h is proportional to the length of the longest edge of the FE mesh.4 .1 Balanced-force implementation As mentioned in the introduction, to avoid dealing explicitly with the discontinuity in Eqs. (6) and (7), aregularized Heaviside function H (cid:15) will be used in practice instead of H in Eq. (5): H (cid:15) ( φ ) = , φ > (cid:15), (cid:18) φ(cid:15) + 1 π sin (cid:18) πφ(cid:15) (cid:19)(cid:19) , | φ | ≤ (cid:15), , φ < − (cid:15). (9)As a result, the transition of fluid properties within a small layer of thickness 2 (cid:15) around the interface canbecome smooth. The parameter (cid:15) ∼ O ( h ) should be as small as possible to have an accurate transition ofmaterial properties while it cannot be smaller than h h for P1 FEs [29]).In general [13, 30, 8], the δ function in Eq. (7) is also replaced by a regularized function δ (cid:15) which is definedby the analytical derivative of H (cid:15) , namely δ (cid:15) = d H (cid:15) d φ . This function may be weighed to deal with large ratiosof mass density [10, 23].However, such formulation leads to oscillations in the case of a static bubble or droplet submitted to nogravity [24]. In this case Eq. (4) has no effect and Eq. (1) reduces to ∇ p = σ s κ s n s δ . Regardless of thenumerical method used to estimate the mean curvature, if the discretization space of ∇ p h does not matchthat of n s ( φ h ) δ ( φ h ), numerical errors will induce spurious non-physical currents [24].To avoid these oscillations, we replace n s ( φ h ) δ ( φ h ) by the numerical gradient of the Heaviside function ∇ H (cid:15) ( φ h ), which leads to the following balanced-force implementation of the CSF model: f s ( φ h ) = σ s κ s ( φ h ) ∇ H (cid:15) ( φ h ) . (10)where H (cid:15) ( φ h ) is defined in Eq. (9) and κ s ( φ h ) is computed using Eq. (8). This is a balanced-force implemen-tation because we are using equal-order discretizations for the pressure, the LS function and the regularizedHeaviside function.In order to compute κ s ( φ h ) directly through Eq. (8), φ h needs to be at least piece-wise quadratic. Alterna-tively, we could use integration by part to eliminate the computation of second derivatives of the LS function[18, 19]. Nevertheless, the choice of a piece-wise quadratic LS function is adopted in this paper to improvemass and energy conservation [29]. As a consequence, p h also needs to be piece-wise quadratic and, eventhough RBVMS stabilization is used, v h needs to be at least piece-wise quadratic because of the inf-supcondition. If we wanted to use a higher-order discretization for φ h and obtain an improved mass conservation[15], we could avoid those restrictions by projecting ∇ H (cid:15) ( φ h ) down to the pressure gradient discretizationspace. However, this would deteriorate the accuracy of the discrete CSF model. Consequently, we adopt aP2/P2/P2 discretization. The stabilization of advection-dominated equations, in particular Navier-Stokes equations in the convection-dominated regime, has been a very active research topic for the last decades [4, 6, 5]. The VMS frameworkproposed in 1998 has given interesting insights on how this stabilization could be achieved by approximatingthose equations simultaneously at two scales ( i.e. , coarse and fine) and projecting the fine scale onto thecoarse scale [6]. The coarse one is simply the one discretized by the mesh, while the subgrid scale should bemodeled separately. The difficulty resides in how the fine scale is modeled and how it is projected to stabilizethe coarse scale. This issue can be addressed through the RBVMS framework which has been proposed forincompressible Navier-Stokes equations [7] and later on extended to the simulation of multiphase flow withsurface tension [18, 8].In the RBVMS framework for multiphase flow [3], the velocity, pressure and LS approximations are written5s v ≈ v h + v (cid:48) p ≈ p h + p (cid:48) φ ≈ φ h + φ (cid:48) (11)where { v h , p h , φ h } are the coarse scale solutions and { v (cid:48) , p (cid:48) , φ (cid:48) } are the fine scale solutions. The discretespace of velocity, pressure and LS test functions { w h , q h , ψ h } is denoted as W h . The semi-discrete weak formof Eqs. (1) and (4) follows [8]:Find { v h , p h , φ h } ∈ V d +2 h , such that, (cid:90) Ω h w h .ρ ( φ h ) (cid:18) ∂ v h ∂t + v h . ∇ v h (cid:19) dΩ+ (cid:90) Ω h ∇ w h : µ ( φ h ) (cid:0) ∇ v h + ∇ T v h (cid:1) dΩ − (cid:90) Ω h p h ∇ . w h dΩ − (cid:90) Ω h w h . f g ( φ h )dΩ + (cid:90) Ω h w h . f s ( φ h )dΩ − (cid:88) K ∈T h (cid:90) K ( ρ ( φ h ) v h . ∇ w h ) . v (cid:48) ( v h , p h , φ h )dΩ+ (cid:88) K ∈T h (cid:90) K ρ ( φ h ) w h . ( ∇ v h . v (cid:48) ( v h , p h , φ h )) dΩ − (cid:88) K ∈T h (cid:90) K ρ ( φ h ) ∇ w h : ( v (cid:48) ( v h , p h , φ h ) ⊗ v (cid:48) ( v h , p h , φ h )) dΩ − (cid:88) K ∈T h (cid:90) K ∇ . w h p (cid:48) ( v h , p h , φ h )dΩ = 0 , ∀ w h ∈ W dh , − (cid:90) Ω h q h ∇ . v h dΩ+ (cid:88) K ∈T h (cid:90) K ∇ q h . v (cid:48) ( v h , p h , φ h )dΩ = 0 , ∀ q h ∈ W h , + (cid:90) Ω h ψ h (cid:18) ∂φ h ∂t + v h . ∇ φ h (cid:19) dΩ − (cid:88) K ∈T h (cid:90) K ∇ ψ h . v h φ (cid:48) ( v h , φ h )dΩ = 0 , ∀ ψ h ∈ W h , (12)where we have indicated explicitly the nonlinear dependence of some terms on the solution. The boundaryintegral term is omitted due to homogeneous Neumann boundary conditions. The fine scale stabilizationterms are integrated element-wise, where T h is the FE mesh which is a set of elements spanning the discretedomain Ω h ≈ Ω. The fine scale stabilization term for the conservation equation makes this formulationinf-sup stable and enables us to use equal-order FEs for velocity and pressure. Readers are referred to Ref.[7] for the details regarding the remaining fine scale terms. In the RBVMS framework, the fine scale solutionsare simply modeled as proportional to the residuals of the strong form coarse scale equations: v (cid:48) ( v h , p h , φ h ) = − τ v ( v h , φ h ) ρ ( φ h ) (cid:18) ρ ( φ h ) (cid:18) ∂ v h ∂t + v h . ∇ v h (cid:19) − ∇ .σ ( v h , p h , φ h ) − f g ( φ h ) + f s ( φ h ) (cid:19) ,p (cid:48) ( v h , p h , φ h ) = − τ p ( v h , φ h ) ρ ( φ h ) ∇ . v h ,φ (cid:48) ( v h , φ h ) = − τ φ ( v h ) (cid:18) ∂φ h ∂t + v h . ∇ φ h (cid:19) , (13)6here σ is the Cauchy stress tensor defined in Eq. (2). Here we have again indicated explicitly the dependenceof the fine scale solutions on the coarse scale solutions as well as that of the following stabilization parameters: τ v ( v h , φ h ) = (cid:32)(cid:18) t (cid:19) + v h . Gv h + C I (cid:18) µ ( φ h ) ρ ( φ h ) (cid:19) G : G (cid:33) − ,τ p ( v h , φ h ) = (tr( G ) τ v ( v h , φ h )) − ,τ φ ( v h ) = (cid:32)(cid:18) t (cid:19) + v h . Gv h (cid:33) − . (14)The constant C I depends on the type and order of the FE, ∆ t is the time step, and G = B T B is the elementmetric tensor, with B the transposed inverse of the gradient of the mapping from reference element to meshelement [7]. Note that those stabilization terms are defined and computed at quadrature points within eachelement, contrary to other stabilization methods where they may be defined as constant element-wise [4, 5].For the time discretization, we use the backward Euler scheme which is unconditionally stable but requiresa Newton-Raphson scheme to deal with the nonlinear terms in Eq. (12), the fine scale solutions in Eq. (13)and the stabilization parameters in Eq. (14). In order to avoid the difficulty of computing all analyticalderivatives, we adopt numerical differentiation. This procedure is conducted locally inside each elementby iteratively adding and then removing a small perturbation δ ND to the node-wise values of each of thesolutions { v h , p h , φ h } . The drawback of this approach is that the Newton-Raphson scheme will not convergeup to a very low tolerance depending on the value of δ ND .If the prescribed time step ∆ t max is too big or the solution is irregular, the Newton-Raphson schememight not converge in a prescribed maximum number of iterations or even diverge and obtain solutionswhich lead to a larger residual from one nonlinear iteration to another. If such situation arises, we restartthe Newton-Raphson scheme with half the time step and repeat this procedure until we find a time step smallenough for the Newton-Raphson scheme to obtain a converged solution in less than 10 iterations. Inversely,when a time increment can be solved at once with no restart, we multiply the time step by 1.5, unless it leadsto a larger value than ∆ t max . In all simulations presented in this paper, this automatic time step controlalgorithm was always successful in making the Newton-Raphson scheme converge in less than 10 iterations.As observed in Sec. 4, this automatic variation of the time step did not lead to an increase of the numberof time increments higher than a factor of two. The use of the regularized Heaviside function as defined in Eq. (9) requires φ to remain a signed distancefunction throughout the simulation. As mentioned in the introduction, the distance property is lost duringinterface motion and deformation, and should be restored using LS reinitialization. Because the implicitscheme enables the use of a large prescribed time step ∆ t max , we reinitialize φ h at every time step.Recently developed indirect methods based on a nonlinear elliptic equation are promising because, in theory,this nonlinear elliptic equation can be easily solved using an FE method with any FE type or order [22].Additionally, it should be theoretically possible to integrate this approach directly into the fully implicitRBVMS framework as an additional regularity constraint [31] and even to add a mass conservation constraintto reduce mass loss during advection and reinitialization [31].In the present work, we adopt a more straightforward approach as we opt for direct geometric reinitialization[20]. The implementation of a P2 LS function is technically challenging because the intersection of eachelement of the mesh with the zero isolevel of φ h is not necessarily planar. To avoid this problem, wehierarchically subdivide each P2 FE twice and reconstruct a fine surface mesh of the zero isolevel of φ h as aP1 mesh of straight line segments (2D) or flat triangles (3D). More details on this method can be found inRef. [23].The distance of any point to the interface Γ , is approximated by looking for the closest projection amongall elements of the surface mesh [32]. This search can be conducted efficiently using a space partitioning7echnique [33, 20]. In a recent work [34], it has been shown that Graphics Processing Units (GPUs) canbe much more efficient for this kind of task than conventional CPU implementations. Therefore, we opt fora brute force search where we project all mesh nodes on all elements of the surface mesh, with no spacepartitioning at all. This search is conducted independently for each mesh node on a separate GPU thread.Additionally, we let the LS reinitialization procedure automatically truncate the LS function around theinterface so that its values be restricted to [ − (cid:15), (cid:15) ]. Thanks to the higher-order discretization of the solution { v h , p h , φ h } , it is expected to improve the accuracyand the conservation of the LS function at each time increment compared to P1 methods. The accuracyshould be controlled near to the interface, in the region of thickness 2 (cid:15) where we operate the fluid propertiestransition defined in Eqs. (9) and (5) and integrate the surface tension force defined in Eq. (10).The challenge is that the interface is evolving during the simulation. We hence propose to dynamically adaptand refine the unstructured mesh during the simulation to track the interface and keep elements refined inthe transition region. To reach this objective, we define a sensor variable s h = ( H (cid:15) ( φ nh ) , H (cid:15) (2 φ nh − φ n − h )),where φ nh is the discrete LS function at the start of a time increment (before Newton-Raphson solution) and2 φ nh − φ n − h is an extrapolation of φ n +1 h . We do not include the velocity and pressure fields in the sensor vari-able in order to avoid difficulties in taking into account errors of different natures. This will be investigatedin a future work as it could be relevant for multiphysics problems. It can be noted that, because of the useof an equal-order P2/P2/P2 FE method, the convergence rate for the velocity and pressure fields should beimproved in the transition region.Isotropic mesh adaption can be based on any scalar a priori or a posteriori estimation of the approxima-tion error at integration points or mesh nodes. For instance, the well-known Zienkiewicz-Zhu error estimator[35] considers the point-wise difference between the numerical solution and the improved approximationrecovered using an operation called Superconvergent Patch Recovery (SPR). This SPR operation which isused subsequently in this work, consists in recovering a higher-order and higher-regularity approximationof the solution around each mesh node. For a current P1 approximation at a given node, the recoveredapproximation should be of order 2 [36], namely it should embed a recovered node-wise continuous gradientand Hessian matrix. This recovered approximation is fitted in a least-squares sense to the values of thenumerical solution at neighboring mesh nodes which is called the patch.Alternatively, the fine scale solutions of the VMS framework described in Subsec. 3.2 can be used as esti-mates of the error in the coarse scale approximations [6, 7].The latter interpolation error can be used as a target for error estimation because the approximation erroris bounded by it. In particular, the interpolation error for a P1 solution can be estimated from the Hessianmatrix of this solution and that for a P2 solution from its third derivatives [37] which can be recoverednode-wise using a SPR one order higher than that of the derivatives to be recovered [38], namely a SPR oforder 3 for second derivatives and order 4 for third derivatives.Elements can be refined and coarsened based on the point-wise estimation of the approximation or interpo-lation error.For anisotropic mesh adaption, additional directional information on the error is required not only torefine and coarsen elements, but also to stretch them differently along each direction. Improved convergencerates are obtained with unstructured anisotropic mesh adaption, because the mesh can be distorted andoriented locally depending on the estimated error [38, 39, 37].In the P1 case, the interpolation error estimator reduces to a Hessian-based error estimator [40, 41, 16]. Aparticular feature of this estimator is that the Hessian matrix defines locally a scalar product and a distancemeasure which are distorted compared to the Euclidean ones [42]. This distortion or stretching depends onthe eigenvalues and eigenvectors of the Hessian matrix. The eigenvectors can also be regarded as the axis ofan ellipsoid and the square root of the inverse of the eigenvalues as the radii of this ellipsoid [43]. An edge-based version of the Hessian-based error estimator has been proposed recently and applied to multiphase8ow simulations [39, 1].To the best of the authors’ knowledge, there has been no anisotropic error estimator based on the fine scalesolutions of the VMS framework in the literature except the anisotropic versions of the Zienkiewicz-Zhuestimator [44, 45]. Anisotropic a priori error estimators for multiphase flow simulations directly targetingthe error in the geometric representation of the interface through the LS function have been proposed inRefs. [46, 47]. In the following, we opt for an error estimator targeting the interpolation error of the sensorvariable. This will allow us not only to adapt the mesh during the simulation but also pre-adapt the meshbefore solving any equation so as to ensure that the interpolation of the initial LS function is accurate.To prepare an unstructured mesh for the initial time increment or any current time increment, it isrequired to estimate the interpolation error of s h , and then adapt the mesh to distribute this error uniformlyon the domain, for a given complexity. Indeed, we will not try to prescribe an error tolerance, which mayblow up the computational cost if the interface area increases during the simulation. Instead, we keep thecomplexity (the number of P1 nodes of the mesh) constant during the simulation and expect the error esti-mator to capture only a certain level of detail which can be afforded with the prescribed complexity. Suchapproaches have been proposed for instance in Refs. [40, 39] for anisotropic P1 interpolation.Approaches to dynamically adapt an unstructured mesh to control the interpolation error for P1 FEs havebeen widely discussed in the literature whereas there has been few works on higher-order discretizations,with the exception of Refs. [37, 48].Recently, the continuous mesh framework for metric-based mesh adaption [49, 50] has been extendedto the higher-order case [37, 48]. In the following, we briefly describe this framework and our proposedimplementation. We denote by s i the component i = 1 . . . dim ( s ) of the continuous sensor variable (dim ( s ) =2 in this work) and by Π s i its discrete P2 interpolation. The approximation error || s i − s h,i || L ( K ) is boundedwithin an element K ∈ T h by the interpolation error || s i − Π s i || L ( K ) which is itself bounded by the L normof the third derivatives ||∇∇∇ s i || L ( K ) : || s i − Π s i || L ( K ) ≤ C ( K ) ||∇∇∇ s i || L ( K ) , ∀ K ∈ T h , ∀ i = 1 . . . dim ( s ) . An approximation of the third derivatives of s i is recovered node-wise using SPR which is naturally extendedfor the P2 sensor variable by recovering the continuous first, second, third and fourth derivatives of s i aroundeach node [38]. Hence, 15 unknowns are solved at each node in 2D and 35 in 3D. To ensure that this problemis well-defined, we use patches with at least twice as many nodes as unknowns for the least-squares fit. Thenode-wise continuous third derivatives of the two-dimensional sensor variable are extracted from this SPRoperation.The main challenge with respect to the P1 case is that ∇∇∇ s i is not a symmetric positive definite second-order tensor and thus does not define a quadratic form, a scalar product and a distance measure. The ideaproposed in Ref. [37] is to find a suitable symmetric positive definite matrix Q ( x ) at each point x ∈ Ω suchthat |∇∇∇ s i ( x ) . y . y . y | ≤ | y T Q ( x ) y | , ∀ y ∈ R d , ∀ i = 1 . . . dim ( s ) . (15)In order to have a close bound, the determinant of Q ( x ) should be as small as possible. This constrainedminimization problem hence consists in finding a symmetric positive definite matrix Q ( x ) such that det Q ( x )is minimal and the bound in Eq. (15) is satisfied. It was solved in Ref. [37] using a so-called logarithmtransformation and a simplex minimization algorithm. Alternatively, a simpler approach consisting in di-rectly finding a quadratic form approximating ∇∇∇ s ( x ) through a least-squares fitting was proposed in Ref.[48]. In this work, we develop a similar approach which consists in directly approximating ∇∇∇ s ( x ) as aquadratic form through the following geometric averaging operation: Q ( x ) = exp s ) d dim ( s ) (cid:88) i =1 d (cid:88) j =1 log (cid:32)(cid:18) ∇∇ ∂s i ∂x j ( x ) (cid:19) − (cid:33) − . (16)9his geometric averaging operation is widely used to average or interpolate metric fields in various applica-tions [42, 43]. It is called geometric because ∇∇ ∂s i ∂x j ( x ) is a symmetric positive definite matrix and can beassociated to an ellipsoid whose radii are given by the eigenvalues of ( ∇∇ ∂s i ∂x j ( x )) − . It is shown in Sec.4 that this simple averaging operation preserves the directions embedded in ∇∇∇ s ( x ) and, in particular,predicts correctly the decrease of the interpolation error in regions with flat interface.The conversion of the directional errors carried by the metric field Q into a metric field M of directional meshsizes requires the solution of a constrained minimization problem consisting in minimizing the total errorwhile uniformly distributing local errors and controlling the complexity. The solution of this constrainedminimization problem in the higher-order case can be expressed as [37] M ( x ) = N d (cid:18)(cid:90) Ω (det( Q ( x )) d d x (cid:19) − d (det( Q ( x ))) − d Q ( x ) , (17)where N is the prescribed number of P1 nodes. It is reminded that the asymptotic convergence rate in L norm with respect to the mesh size for a uniform isotropic P2 mesh is 3, which means a rate of 3 d withrespect to N .The metric field defined in Eq. (17) can be used to compute the distorted lengths of each edge of thecurrent FE mesh and determine whether they are too long (distorted length larger than one) or too short(distorted length smaller than one). The distorted volumes of each element of the current mesh can alsobe computed. An optimal mesh for a metric field M has a minimum number of elements whose distortedvolumes are as uniform as possible and edges whose distorted lengths are as close as possible to one.To solve this multi-objective optimization problem, we use a mesh adaption algorithm developed in previousworks [51, 9]. Note that we keep the P2 nodes at edge centers in this optimization and the unstructuredmesh remains linear (only the interpolation is quadratic). For a given element, the quality criterion which isused in the algorithm considers the element’s distorted volume and the average of the distorted lengths of theelement’s edges. All conformal topological modifications including node re-positioning and edge splitting aretested in small patches in the neighborhood of each node and edge of the mesh, and the modification whichmaximizes the quality criterion is accepted. This iterative algorithm with local modifications ends when nomore local modification can be found to improve the quality of the mesh. The details on the mathematicalbackground of this algorithm can be found in Ref. [51], while the implementation details can be found inRef. [9].Once the mesh adaption algorithm terminates, the velocity field and the LS function should be transferredfrom the old mesh to the new (adapted) one. First, a space partitioning technique is used to locate efficientlythe element of the old mesh containing each node of the new mesh. The values of the velocity and the LSfunction are then computed at each node of the new mesh using P2 interpolation from the containing elementof the old mesh. At this step, the P2 discretization of both the velocity field and the LS function can reduceany errors and diffusion induced by mesh adaption and fields transfer. Finally, the new mesh and thetransferred velocity field and LS function can be used for the current time increment. The proposed method has been implemented in a C code (as of Ref. [52]), whereas the LS reinitializationprocedure mentioned in Subsec. 3.3 has been implemented in NVIDIA’s CUDA programming language forGPU computing and the mesh adaption algorithm presented in Subsec. 3.4 has been implemented in
C++ (as of Ref. [53]). For computationally demanding operations such as numerical differentiation and integrationof the discrete RBVMS weak form presented in Subsec. 3.2, implementations rely on shared-memory multi-threaded parallelism using
OpenMP [54]. The lis library to solve the linear system at each iteration of theNewton-Raphson scheme also relies on
OpenMP [55]. In particular, we employ the generalized minimum10esidual method preconditioned using an incomplete LU factorization with threshold. The tolerance forthis linear solver is set to 10 − , the perturbation for numerical differentiation δ ND is set to 10 − , and thetolerance for the Newton-Raphson solver is set to 10 − .Following the notations in Ref. [56], we note NDOF the number of degrees of freedom for a given simulation.It is reminded that depending on the boundary conditions, there are at most d + 2 DOF per node of the FEmesh. The number of time increments for a simulation is denoted as NTS, and the number of elements ofthe mesh as NEL. This first set of simulations is conducted to assess the efficiency of the balanced-force implementation of theCSF model as well as the accuracy of curvature computation using P2 FEs in comparison with other works[57, 24, 8, 26]. Following these works, we define the computation domain as a cube of dimensions 8 × × × ×
40 and 20 × ×
20 elements. A spherical droplet of radius R = 2 isplaced at the center of the computation domain as shown in Fig. 1(a) and no-slip conditions are applied atall the computation domain boundaries. Non-dimensionalized fluid and flow properties ρ = 1, ρ = 10 − , µ = µ = 0, g = 0, σ s = 73 are used. The simulation is conducted for 50 time increments using a time step∆ t = 0 . p between outside and inside the droplet should converge to that given by the Young-Laplace equation:∆ p = σ s κ s = 73 , with κ s = 2 R = 1 . As an error measure, we use the maximum Euclidean norm of the velocity vector in the domain max x ∈ Ω ( || v ( x , t ) || ),after one time increment and after 50. At the same instants, we also report the error of the difference betweenthe highest and the lowest pressures Error(∆ p max ( t )). (a) (b) Figure 1: 3D single static inviscid droplet: (a) uniform 40 mesh with droplet placed at its center, (b)pressure field computed at t = 0 .
050 using the 20 uniform P2 mesh and exact curvature.First, using the exact curvature in the numerical scheme, both errors should be below the tolerance usedfor the Newton-Raphson scheme, owing to the balanced-force implementation of the CSF model. In this11ase, we can also conduct the simulation using equal-order P1 FEs, because we do not need to computethe curvature. After 50 time increments, errors might increase depending on the accuracy of the numericalmethod, in particular the LS reinitialization method.Results are reported in the first part of Tab. 1. At equal NDOF, both errors of pressure and of velocityare below the tolerance used for the Newton-Raphson scheme, even after 50 time increments, regardless ofthe order of the FE method. This result is in agreement with the theory and validates our balanced-forceimplementation. Contrary to other works [26], we cannot reduce this error down to machine precision, whichis mainly due to the numerical differentiation that uses δ ND = 10 − .Note that in Tab. 1 we have also reported the regularization thickness (cid:15) . There can be significant oscillationsof the pressure field near the interface if (cid:15) is much smaller than the mesh size. For instance, the pressurefield obtained using exact curvature and a value of (cid:15) corresponding to a regularization over four elements isshown in Fig. 1(b). Mesh 40 (cid:15) . . . . p max ( t = 0 . < − < − . × − . × Error(∆ p max ( t = 0 . . × − . × − . × . × max x ∈ Ω ( || v ( x , t = 0 . || ) < − < − . × − . × − max x ∈ Ω ( || v ( x , t = 0 . || ) 8 . × − . × − . × . × Table 1: 3D single static inviscid droplet simulation results for different FE meshes and orders, using eitherexact or numerical curvature. Errors below δ ND = 10 − are not reported.Second, using numerically computed curvature, errors will depend on the interpolation error of the LSfunction and its second derivatives. In the literature, most authors rely on derivatives recovery techniqueswhich can smooth the mean curvature [57, 24, 25, 26, 3]. A higher-order and higher-regularity isogeometricanalysis method has also been proposed [8]. Here, we use a P2 discretization of the LS function withoutderivatives recovery nor smoothing technique.Results are reported in the second part of Tab. 1. There are clearly numerical errors in the computation ofthe mean curvature, which may significantly affect simulation results depending on the discretization at thevicinity of the interface. These errors accumulate with time steps and are hence higher after 50 steps. Thiserror accumulation effect seems to be alleviated with a second order time stepping scheme [8] which couldbe worth considering.In the literature [57], various derivatives recovery and smoothing techniques using the finite difference methodhave been reported to lead to errors of the velocity after one time increment between 8 . × − and3 . × − for the 40 mesh. After 50 time increments, those errors have been reported to vary between3 . × − and 2 . × depending on the smoothing technique. These works did not employ a balanced-force algorithm. For example, the velocity errors between 4 . × − and 4 . × − after one timeincrement, and between 4 . × − and 1 . × − after 50 time increments have been reported in Ref.[24] using the finite difference method with different smoothing techniques for the 40 mesh. Overall, ourresults after one time increment are better than those obtained with non balanced-force implementations,whereas they are not as accurate as those reported in Ref. [24]. It seems that the balanced-force implemen-tation is necessary even though it is not sufficient to eliminate spurious currents. Thus, the balanced-forceimplementation should be combined with curvature smoothing techniques.This speculation is confirmed by the fact that our errors are at least two orders superior to those reportedin Ref. [8] using equal-order P1 FEs and a projection technique to recover the mean curvature. Althoughsuch technique would make the numerical differentiation in our fully implicit scheme non-local and hence12ore difficult to implement, it may be worth considering in the future.Ref. [8] reported that isogeometric analysis with its higher-order basis functions led to a more accurate meancurvature computation but did not lead to a reduction of spurious currents, compared to a P1 discretization.This is also confirmed by the results using numerical curvature in Tab. 1.However, higher-order methods are still relevant to improve the accuracy on the pressure field. Errors after50 time increments of 1 . × for the 40 mesh and 3 . × for the 20 mesh have been reported in Ref.[8] using the equal-order P1 FE method with projection, and 7 . × − and 1 . × using isogeometricanalysis. Our results are of the same order but we observe a drastic increase of the error after 50 timeincrements, which emphasizes again the need to improve the time stepping scheme in future developments.As a conclusion, this first set of simulations validates our balanced-force implementation of the CSFmodel, because all computations using exact mean curvature lead to errors both in velocity and pressurewhich are below the tolerance used for the Newton-Raphson solver. Regarding computations with numericallycomputed mean curvature, spurious non-physical currents are not completely eliminated. However, theproposed equal-order P2 framework enables a good compromise between the accuracy of the pressure fieldand the simplicity of the scheme because we do not use any projection or curvature smoothing technique. This second set of simulations considers the dynamics of bubble rise in a 2D domain. The simulation setupshown in Fig. 2 has been proposed in a benchmark for multiphase flow simulation methods [56]. Thecomputation domain is a rectangle of dimensions 1 ×
2, discretized with meshes of 40 ×
80, 80 ×
160 and160 × R = 0 .
25 is initially placed at (0 . , . • for case 1, properties are ρ = 1000, ρ = 100, µ = 10, µ = 1, g = 0 . σ s = 24 . • and for case 2, ρ = 1000, ρ = 1, µ = 10, µ = 0 . g = 0 . σ s = 1 . L errors of thebubble area V b , the bubble center of mass y coordinate y b and the bubble rise velocity v b integrated overtime. Those parameters are respectively given by V b ( t ) = (cid:90) Ω H ( − φ ( x , t ))dΩ ,y b ( t ) = (cid:82) Ω H ( − φ ( x , t )) x dΩ V b ( t ) . e y ,v b ( t ) = (cid:82) Ω H ( − φ ( x , t )) v ( t )dΩ V b ( t ) . e y . The relative L error integrated over time isError( q ) = (cid:118)(cid:117)(cid:117)(cid:116) (cid:82) T ( q ( t ) − q ref ( t )) d t (cid:82) T ( q ref ( t )) d t , where T = 3 and q = V b , y b , or v b . Simulations with different space (NDOF) and time (NTS) discretizationsare conducted with and without anisotropic mesh adaption. The reference solution V b,ref is πR , as theflow is incompressible. Since we do not have analytical solutions for y b,ref and v b,ref , we use the numerical13olutions corresponding to the simulation with best V b accuracy for each case.The asymptotic convergence rate with respect to the square root of NDOF is 3, as explained in Subsec.3.4, while the asymptotic convergence rate with respect to NTS is 1. Even though the numerical schemeis unconditionally stable, these optimal convergence rates will not be reached if the mesh size is too smallcompared to the time step. In Ref. [56], convergence rates corresponded to a simultaneous reduction of themesh size and the time step. Furthermore, for a given NDOF, increasing NTS may not necessarily increasethe accuracy as more time increments mean more LS reinitializations. Indeed, the LS reinitialization method(Subsec. 3.3) slightly shifts the interface at each reinitialization. g Fluid 2 (gas)Fluid 1 (liquid)0.5 0.5 R=0.25 0.5 2 v x =v y = v x =v y = v x = v x = Figure 2: Setup for the 2D single bubble rise simulations.
In case 1, the bubble which was initially circular becomes nearly elliptical with a flat bottom surface. One ofour solutions computed using an adaptive mesh is shown in Fig. 3. Anisotropic mesh adaption is particularlyrelevant for this kind of interface morphology because elements can be stretched along directions where thethird derivatives of the LS function are close to zero, to reduce the computational cost. This can be seen inFig. 3(c) in the lower part of the bubble. This anisotropic stretching of elements is also relevant to controlthe computational cost, because there is a longer perimeter to cover in Fig. 3(c) than in Fig. 3(a), while thenumber of elements is almost the same, as explained in Subsec. 3.4.A closer look at the final shapes of the bubble for some of the simulations with adaptive mesh is presentedin Fig. 4. A good convergence is obtained by dint of the fine mesh which is maintained near the interface,using anisotropic mesh adaption.Errors in the computation of the parameters of interest are reported in Tab. 2. Although NTS wasprescribed to be equal to 100, 200, 400 and 800, the Newton-Raphson problem was too stiff for the simulations14 a) (b) (c)
Figure 3: Velocity field in rise direction v . e y for the 2D single bubble rise case 1 simulation using an adaptivemesh with NDOF ≈ t = 0, (b) t = 1 .
5, (c) t = 3 . NDOF ≈ ≈ ≈ Figure 4: Liquid-gas interface for some of the 2D single bubble rise case 1 simulations with adaptive meshat t = 3. 15ith very fine meshes for the given time step. Thus, the time step was automatically decreased to obtain asolution as explained in Subsec. 3.2. The most extreme case was the simulation with uniform 80 ×
160 meshand a prescribed NTS of 100 which was automatically increased to nearly 200. Results from this simulationwere not reported in Tab. 2 to avoid any confusion.Mesh 20 ×
40 40 ×
80 80 ×
160 Adaptive Adaptive Adaptive (cid:15) . .
05 0 .
025 0 .
05 0 .
025 0 . ≈ ≈ ≈ ≈ ≈ ≈ V b ) 5 . × − . × − . × − . × − . × − Error( y b ) 4 . × − . × − . × − . × − . × − Error( v b ) 8 . × − . × − . × − . × − . × − NTS 207 200 200 200 204 261Error( V b ) 3 . × − . × − . × − . × − . × − . × − Error( y b ) 4 . × − . × − . × − . × − . × − . × − Error( v b ) 7 . × − . × − . × − . × − . × − . × − NTS 406 400 400 400 407 463Error( V b ) 2 . × − . × − . × − . × − . × − . × − Error( y b ) 7 . × − . × − . × − . × − . × − . × − Error( v b ) 8 . × − . × − . × − . × − . × − . × − NTS 802 800 800 800 804 880Error( V b ) 3 . × − . × − . × − . × − . × − . × − Error( y b ) 1 . × − . × − . × − . × − . × − Error( v b ) 7 . × − . × − . × − . × − . × − Table 2: Relative L errors for 2D single bubble rise case 1 simulations using different space (NDOF) andtime (NTS) discretizations, with and without anisotropic mesh adaption. All results with convergence ratewith respect to the square root of NDOF higher than 2.5 are highlighted in bold.As shown in the left part of Tab. 2 where no mesh adaption was used, there is a compromise to findbetween the time discretization and the space discretization. Optimal convergence rates in V b and y b arereached for instance for the intermediate 40 ×
80 mesh using the finest time discretization (highlighted inbold in Tab. 2). This confirms that if the mesh is not fine enough, interface shifts due to LS reinitializationbecome the dominant source of error and lead to an increased error with smaller time steps.Overall good mass conservation is obtained because of the use of a P2 LS function. This is drasticallyimproved using an adaptive mesh. For instance, using an adaptive mesh with NDOF ≈ ,
000 (right partof Tab. 2), the error of V b for NTS = 400 is lower than the errors obtained using a uniform mesh withsame NTS, even using the finest uniform mesh. A factor of reduction in terms of NDOF close to 35 is henceachieved.Anisotropic mesh adaption is also relevant to reach optimal convergence rates for a lower NDOF than uni-form meshes. This can be seen in the last line (NTS ≈ V b is obtained using an adaptive mesh with NDOF ≈ , ≈ , .2.2 Case 2 In case 2, thin bubble filaments develop at lateral extremities of the bubble and distort significantly, whilethe main part of the bubble remains nearly elliptical with a flat bottom surface. One of the solutions ob-tained using an adaptive mesh is shown in Fig. 5. The great advantage of anisotropic mesh adaption isdemonstrated once again, because both thin filaments are captured and tracked by the adaptive mesh, asshown in Fig. 5(f).A closer look at the final shapes of the bubble in some simulations is presented in Fig. 6. Convergenceis more difficult for this case as the filaments become thinner and thinner with a smaller mesh size and asmaller (cid:15) . This change of thickness of the bubble filaments with finer meshes affects slightly the whole shapeof the bubble. For the simulation with the highest NDOF and NTS, the lateral extremities of the bubbleseem to transform into new bubbles which could potentially emerge and breakup. (a) (b) (c)(d) (e) (f)
Figure 5: Velocity field in rise direction v . e y for the 2D single bubble rise case 2 simulation using an adaptivemesh with NDOF ≈ t = 0, (b) t = 1 .
5, (c) t = 3 .
0. Zoomed images on theadapted mesh at: (d) t = 0, (e) t = 1 .
5, (f) t = 3 .
0. 17
DOF ≈ ≈ ≈ ≈ Figure 6: Liquid-gas interface for some of the 2D single bubble rise case 2 simulations with adaptive meshat t = 3. 18ome codes assessed in the benchmark predicted a breakup of the thin filaments [56]. This could beclearly observed by looking at an additional parameter of interest, viz. the bubble circularity S b which isdefined as S b ( t ) = 2 πV b ( t ) (cid:115) (cid:82) Ω δ ( − φ ( x , t ))dΩ π . After the point of breakup, some codes did not show clear convergence on circularity curves and significantdifferences were found when comparing curves computed using different codes. We provide an explanationof those differences in the following.The TP2D code used the LS method with bilinear FEs for the LS function and a second order time step-ping scheme with a weak coupling between Navier-Stokes equations and the LS advection equation. TheFreeLIFE code also used the LS method with a second order time stepping scheme and a weak coupling,but with linear FEs for the LS function. An additional algorithm was implemented to ensure global massconservation. The MooNMD code used linear FEs with conformal meshing of the interface and a secondorder time stepping scheme with a weak coupling. An Arbitrary Lagrangian-Eulerian (ALE) method wasimplemented to achieve local mass conservation better. See Ref. [56] for more references on each code andhow they were used in the benchmark.Based on these facts, the code that could theoretically achieve the best performance regarding local massconservation was the MooNMD code owing to its ALE method. It is also the only code that could reach aclear convergence on circularity curves for the numerical parameters used in the benchmark. Our methodcan achieve a similar performance by virtue of the use of a P2 LS function and anisotropic mesh adaption.This is demonstrated by looking at the errors of V b in Tab. 3 and for case 1 in Tab. 2.Neither the MooNMD code nor our method predict breakup of the filaments. Moreover, as shown in Fig.7(a), the circularity evolution predicted by our method converges at t = 3 to a value a little greater than 0 . .
2, which matches with the results found in the benchmarkfor all codes. The convergence is nevertheless not as fast as for the circularity, as shown in Tab. 3 where theresults with optimal convergence rate are highlighted in bold. We do not obtain any convergence rate onthe square root of NDOF which is close to three for the rise velocity. A multi-criterion error estimator foranisotropic P2 mesh adaption which would take into account the velocity field should hence be consideredin future developments.We note however that we used NTSs which were lower than that used by the MooNMD code (between 3,000and 6,000). This is a limitation of the ALE approach, which is avoided with the Eulerian method usedherein with anisotropic mesh adaption. Owing to the P2/P2/P2 scheme, we can reach a high NDOF witha relatively low NEL compared to that used by all codes assessed in the benchmark. This is of particularinterest for our approach regarding mesh adaption, because the computational cost of mesh adaption isrelated to NEL (Subsec. 3.4).Following the diagonal in Tab. 3, the total computation times on a machine with a 20-core
Intel XeonE5-2630 v4 @ 2.20GHz
CPU and a 3072-core
NVIDIA GTX TITAN X
GPU were about 10min, 28min,1h28min and 4h32min, respectively. For all case 1 and case 2 simulations, the computation time was mainlyspent on the Newton-Raphson solution of the RBVMS formulation whereas the computational cost of LSreinitialization was negligible (smaller than 1%). For the simulations with an adaptive mesh, the computationtime of error estimation, SPR, metric computation, mesh adaption and fields transfer represented between1 and 15% of the total computation time. In addition, we did not observe any significant mass loss due toP2 fields transfer after mesh adaption. Instead, we observed an improved mass conservation at equal NDOFwith P2 mesh adaption. Therefore, we can conclude that there is a great merit in using an adaptive P2mesh. 19 .
05 0 .
025 0 . . ≈ ≈ ≈ ≈ ≈ ≈ ≈ ≈ V b ) 4 . × − . × − Error( S b ) 2 . × − . × − Error( y b ) 1 . × − . × − Error( v b ) 8 . × − . × − NTS 203 236 323Error( V b ) 1 . × − . × − . × − Error( S b ) 2 . × − . × − . × − Error( y b ) 1 . × − . × − . × − Error( v b ) 8 . × − . × − . × − NTS 403 419 530Error( V b ) 3 . × − . × − . × − Error( S b ) 2 . × − . × − . × − Error( y b ) 1 . × − . × − . × − Error( v b ) 8 . × − . × − . × − NTS 804 825 906 1447Error( V b ) 1 . × − . × − . × − . × − Error( S b ) 2 . × − . × − . × − Error( y b ) 1 . × − . × − . × − Error( v b ) 7 . × − . × − . × − Table 3: Relative L errors for 2D single bubble rise case 2 simulations using different space (NDOF) andtime (NTS) discretizations, with anisotropic mesh adaption. All results with convergence rate with respectto the square root of NDOF higher than 2.5 are highlighted in bold. (a) (b) Figure 7: Circularity (a) and rise velocity (b) evolution for the 2D single bubble rise case 2 simulations usingan adaptive mesh with various NDOF and NTS. 20 .3 Single bubble rise (3D)
Different 3D versions of the single bubble rise benchmark have been proposed in the literature [61, 62, 63, 26].Here we opt for a simple 3D extension where the computation domain is a box of dimensions 1 × × R = 0 .
25 is initially placed at (0 . , . , . ≈ , ≈ ,
000 for a NTS ≈ ≈ ,
000 andNDOF ≈ ,
000 for a NTS ≈ (cid:15) is 0 .
05 for the coarse discretization and 0 .
025 for the fineone.The results with the fine discretization are shown for case 1 in Fig. 8 and for case 2 in Fig. 9. The final shapeof the bubble for each case is similar to that obtained by other authors in the literature [61, 26]. Contraryto the 2D version, there is no excessive stretching of the bubbles into filaments that might breakup for case2, which is also in agreement with the literature [61, 26].A closer look at the final shapes is shown in Fig. 10. There is clearly an overestimation of the bubblevolume in both cases for the coarse discretization, which can be attributed to the low order of the time step-ping scheme based on the 2D simulations. The final shapes are nevertheless smooth, even using the coarsediscretization. In comparison with other results for a similar number of elements found in the literature[61, 26, 8], this result shows the great advantage of using an anisotropic adaptive P2 interpolation for theLS function. (a) (b) (c)
Figure 8: Velocity field in rise direction v . e y for the 3D single bubble rise case 1 simulations using an adaptivemesh with NDOF ≈ ≈
200 at: (a) t = 0, (b) t = 1 .
5, (c) t = 3 . V b ( t ) now represents21 a) (b) (c) Figure 9: Velocity field in rise direction for the 3D single bubble rise case 2 simulations using an adaptivemesh with NDOF ≈ ≈
200 at: (a) t = 0, (b) t = 1 .
5, (c) t = 3 . S b is replaced by the bubble sphericity given by: S b ( t ) = 4 πV b ( t ) (cid:32) (cid:115) π (cid:90) Ω δ ( − φ ( x , t ))dΩ (cid:33) . Time evolution of each parameter of interest is reported for both cases in Fig. 11. For case 1, the errors of V b are 1 . × − for the coarse discretization and 5 . × − for the fine discretization, while for case 2they are respectively 1 . × − and 5 . × − . As expected, the errors are slightly higher for case 2, whilethis difference is quite low compared to the differences between the two simulation results shown in Fig. 10.This proves that our numerical method is robust and can achieve a similar accuracy for very different bubbledeformations.The error of y b is very low for both cases, as shown in Fig. 11(c,d) where the overall position of the bubbleis well captured. In order to describe the bubble’s shape, a finer discretization is required, as shown by thesphericity curves in Fig. 11(a,b). Small oscillations can be seen for case 1 in Fig. 11(a), which are due tolocal mass loss during fields transfer after mesh adaption. Similar oscillations can be seen for case 2 in Fig.11(b). Their amplitude is clearly low, however compared to the overall evolution of the sphericity. The risevelocity is underestimated with the coarse discretization, as shown in Fig. 11(e,f). This was also observed inRef. [26]. The differences between the two discretizations are however quite low. Both the instant at whichthe maximum rise velocity is reached and the value of the maximum rise velocity are well predicted with thecoarse discretization.All 3D single bubble rise simulations were run on a machine with a 20-core Intel Xeon E5-2630 v4 @2.20GHz
CPU, and a 3072-core
NVIDIA GTX TITAN X
GPU. The total computation time was about 9h forcoarse discretization simulations and 34h for fine discretization simulations. About 65% of this computation22 a) (b) (c)(d) (e) (f)
Figure 10: Shapes with transparent interface at t = 3 for the 3D single bubble rise simulations: (a) case 1using the coarse discretization, (b) case 1 using the fine discretization, (d) case 2 using the coarse discretiza-tion, (e) case 2 using the fine discretization. Midsection of the interface using the coarse (orange) and fine(green) discretizations: (c) superposition of (a) and (b) for case 1, (f) superposition of (d) and (e) for case 2.23 a) (b)(c) (d)(e) (f) Figure 11: 3D single bubble rise simulation results using an adaptive mesh with various NDOF and NTS.Case 1: (a) sphericity, (c) center of mass y coordinate, (e) rise velocity. Case 2: (b) sphericity, (d) center ofmass y coordinate, (f) rise velocity. Scales may differ from left to right.24ime was spent on Newton-Raphson solution, 7% on the mesh adaption algorithm, 23% on fields transferafter mesh adaption and 5% on LS reinitialization. The computation time of the other operations such aserror estimation and the SPR operation was negligible.The computational cost caused by the modification of the mesh at each time increment was clearly moresignificant for these 3D simulations. However, this was mainly due to the fields transfer operation and inparticular to the search of the containing element of the old mesh for each node of the new mesh. This couldbe significantly reduced by the implementation of a parallel version of this search on the GPU, as we did forLS reinitialization.The increase in the computational cost by anisotropic P2 mesh adaption is marginal and acceptable giventhe drastic improvements of accuracy in the 2D simulations. Incompressible two phase flows with surface tension were modeled through the Navier Stokes equations, theLevel-Set (LS) method and the Continuum Surface Force (CSF) model. The equal-order piece-wise quadratic(P2) finite element discretization was used for the simultaneous solution of the velocity, the pressure and theLS function. The Navier-Stokes equations with the CSF model and the LS advection equation were stronglycoupled in a residual-based variational multiscale formulation and solved using a Newton-Raphson schemewith numerical differentiation. The fully implicit formulation was shown to be robust and stable even forlarge time steps.A balanced-force implementation was used for the CSF model. This implementation could be achieved in astraightforward manner with no need to rely on the numerically complex projections of pressure gradient,because the same discretization was used for the pressure field and the regularized Heaviside function. Thisimplementation was shown to reduce errors of both the velocity and pressure fields below the tolerance usedfor the Newton-Raphson solver on a static problem using exact curvature.The original approach based on anisotropic P2 mesh adaption was proposed to automatically refine, stretchand orient elements along the interface. As the interface was represented by the zero isolevel of an LSfunction, the new error estimator targeting the P2 interpolation error of the regularized Heaviside functionwas developed to define a directional error metric field. The continuous mesh framework was used to convertlocal directional errors into local directional mesh sizes with a control of the global complexity. This strategywas shown to be more accurate than a uniform fixed mesh for a same number of degrees of freedom, andto enable reductions of the number of degrees of freedom with a factor up to 15 for a same accuracy. Inaddition, the convergence rates and the mass conservation properties could be improved by virtue of P2interpolation.
Acknowledgments
The authors would like to acknowledge the financial support of the French Agence Nationale de la Rechercheto COMP3DRE project (grant number: Projet-ANR-16-CE08-0042).