A High Order Flux Reconstruction Interface Tracking Method Using Preconditioned Phase Field
AA High Order Flux Reconstruction Interface Tracking Method UsingPreconditioned Phase Field
Jabir Al-Salami , Mohamed M. Kamra , and Changhong Hu ∗ Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, Fukuoka, Kasuga, 816-8580 Japan Research Institute for Applied Mechanics, Fukuoka, Kasuga, 816-8580 Japan
Abstract
This paper presents a simple and highly accurate method for capturing sharp interfaces moving in divergence-free velocity fields using the high-order Flux Reconstruction approach on unstructured grids. A well-knownlimitation of high-order methods is their susceptibility to the Gibbs phenomenon; the appearance of spuriousoscillations in the vicinity of discontinuities and steep gradients makes it difficult to accurately resolve shocks orsharp interfaces. In order to address this issue in the context of sharp interface capturing, a novel, preconditionedand localized phase field method is developed in this work. The numerical accuracy of interface normal vectors isimproved by utilizing a preconditioning procedure based on the level set method with localized artificial viscositystabilization. The developed method was implemented in the framework of the multi-platform Flux Reconstruc-tion open-source code PyFR [1]. Numerical tests in 2D and 3D conducted on different mesh types showed that thepreconditioning procedure significantly improves accuracy. The results demonstrate the conservativeness of theproposed method and its ability to capture highly distorted interfaces with superior accuracy when compared toconventional and high-order VOF and level set methods. The high accuracy and locality of the proposed methodoffer a promising route to carrying out massively-parallel, high accuracy simulations of multi-phase, incompress-ible phenomena.
Keywords—
High-order methods, Interface capturing, Flux Reconstruction, PyFR
Simulating multi-phase phenomena requires an accurate computational fluid dynamics (CFD) scheme that is capable ofcapturing or tracking the interface between the different phases or immiscible fluids. Approaches to this problem broadlyfall within two frameworks:Lagrangian interface tracking methods and Eulerian interface capturing methods . Lagrangianapproaches represent the interface explicitly using either tracking markers or particles (eg. MAC methods [2]), or a continuouslydeforming mesh conforming to the interface (eg. Arbitrary Lagrangian-Eulerian (ALE) methods [3, 4, 5]). While interfacetracking methods are suitable for moderately deforming free-surfaces and fluid-structure interactions, they often struggle withseverely deforming and fragmenting interfaces [6]. Furthermore, a major draw-back associated with such methods is thecomputational cost incurred by frequent re-meshing. In Eulerian methods, the interface is implicitly captured by means of anauxiliary field whose value corresponds to different fluids, and that is transported and deformed by the velocity variation in thedomain. Among such methods, the Volume of Fluid (VOF) and level-set methods are the most successful and widely adopted,having been proven effective in a variety of disciplines such as flooding and violent impact problems [7, 8, 9], combustion[10, 11], atomization and evaporation [12, 13].In VOF methods, the auxiliary variable, α , is the volume fraction. It ranges between zero and one and represents the ratioof primary to secondary fluid volume in a discrete computational cell. Therefore, the interface between two immiscible fluidsis represented by a jump in α . If the volume fraction is advected accurately and the interface thickness is kept constant inspace and time, VOF schemes have been shown to demonstrate a high level of accuracy and mass conservation. In order toachieve that, however, the most accurate VOF schemes rely on complicated and computationally intensive geometric interfacereconstruction operations. Another VOF approach is the tangent hyperbolic interface capturing (THINC). THINC methods ∗ Corresponding author. Email: [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] J un tilize a tangent hyperbolic function to reconstruct the volume fraction field. This allows for controlling and maintaining theinterface thickness using a single parameter [14]. A common drawback of VOF methods is the difficulty of obtaining accurateestimates of interface-normal vectors and interface curvature, which are especially important in surface tension dominatedflows and applications that require setting special boundary conditions at the free surface.The phase field approach to interface tracking [15, 16, 17] is similar to VOF in that the interface is implicitly representedby the sharp transition in the auxiliary variable (the phase field in this case). However, instead of resorting to explicitreconstruction of the interface, a diffusion and anti-diffusion term in the transport equation of the phase field maintain adiffuse interface with controllable thickness.Level set methods utilize a signed distance function, whereby the interface is implicitly defined as the zero-level set contour.The smooth variation of the distance function across the interface allows for accurate curvature and normal vector estimates.A major caveat of this approach, however, is the mass loss due to numerical dispersion. This issue has prompted a greatdeal of research efforts which attempted to enhance mass conservation by a variety of approaches, such as coupling to a VOFmethod [18, 19], casting the level set equations in a conservative form [15, 20], and adaptive mesh refinement (AMR) [21], toname a few.Recent trends in CFD research indicate a steadily increasing interest in high-order numerical methods. The ability of suchmethods to produce results with more accuracy on coarser grids when compared to conventional low order methods resultedin a growing consensus among CFD practitioners that high-order methods may constitute the basis of next-generation CFDresearch tools [22, 23]. In line with this vision, some efforts in multi-phase CFD research aimed at using the superior numericalproperties of high-order methods to address some of the limitations of level set and VOF approaches. Qian et. al used thehigh-order Weighted Essentially Non-Oscillatory(WENO) approach in their coupled THINC/level set scheme to obtain high-order polynomial representations of interfaces [19]. Matsushita and Aoki [24] employed the conservative phase field methodwith a 5th order WENO discretization to carry out weakly compressible free-surface simulations.Attempts to utilize polynomial basis high-order methods primarily focused on leveraging the high accuracy and low numericaldissipation of the discontinuous Galerkin (DG) method to improve level-set results. Jibben and Herrmann [25] developeda Runge–Kutta DG approach to re-initialization by solving banded conservative level sets on a separate, refined Cartesiangrid. This approach, however, produces a smeared interface where shearing or severe deformation occurs. Zhang and Yue [26]developed an AMR-DG method where in interface cells the gradient of the level-set function is determined by a weighted localprojection scheme, but mass-loss and interface fragmentation artifacts are still present. Karakus et al. [21, 27] developed aGPU-accelerated AMR re-initialization approach to level set with artificial diffusion stabilization, yet mass loss remained anissue. To solve the issue of oscillatory level-set in reinitialization, Gross et al. [28] used a combination of artificial diffusionand exponential filtering. Most previous DG interface-capturing research efforts shied away from VOF-like sharp interfacerepresentations. This may be attributed to high-frequency oscillations that arise in the vicinity of sharp jumps due to theGibbs phenomenon, an issue that is emblematic of polynomial-basis high-order methods [29].This paper addresses the aforementioned challenges and proposes a novel high-order method for interface capturing usingthe flux reconstruction (FR) approach. The FR approach was proposed by Hyunh [30], and it encompasses the DG and spectraldifference (SD) high-order approaches. The compact spatial stencil of the FR discretization in addition to its amenability toexplicit time marching makes it particularly suitable for modern high-performance-computing hardware, which is characterizedby an excess of computational power relative to memory bandwidth [1, 31].In order to accurately capture sharp interfaces without introducing numerical oscillations, this work proposes a novelapproach based on the conservative phase field method [17]. The proposed approach leverages the accuracy of the FRdiscretization to capture artifact-free, sharp interfaces while maintaining good mass conservation. In order to alleviate theissue of erroneous interface curvature and normal vectors, this paper proposes a level set-based, nonlinear preconditioningequation with localized artificial viscosity stabilization. This results in a simple to implement algorithm, with minimalcoupling overhead and high level of accuracy. This method, henceforth referred to as the Flux Reconstruction PreconditionedPhase Field (FR-PCPF) method, has been implemented in the framework of the open-source code PyFR [1], which is ahigh-performance, multi-platform flux reconstruction code.This paper begins with a brief review of the level set and phase field methods. This is followed by a description ofspatial and temporal discretization in the framework of the flux reconstruction method and other implementation details. Theproposed method is validated and aspects of its performance are assessed by carrying out 2 and 3-dimensional interfacingcapturing benchmarks. The results are also compared to other previous published studies using low and high-order methods.Finally, concluding remarks are given with a brief overview of future work. Phase field based interface capturing techniques are mostly derived from either the Allen-Cahn [32] equation or Cahn-Hillardequation [33], which were proposed to describe the spontaneous separation in multi-component fluid mixtures. The Allen-Cahnequation is more commonly used to model multi-phase flows since it is simpler to implement and, unlike the Cahn-Hillard quation, it does not contain high-order derivatives. The conventional Allen-Cahn equation reads ∂φ∂t = γ (cid:20) ∇ φ − F (cid:48) ( φ ) (cid:15) (cid:21) (1)where φ is the phase field variable in the range of zero to one, γ and (cid:15) are parameters controlling the width of the interface,and F ( φ ) = (cid:0) φ (cid:0) − φ (cid:1)(cid:1) is the so-called double-well potential. The fact that original equation does not guarantee massconservation motivated attempts to enforce conservation through a variety of approaches; such as using Lagrange multipliers[34, 35] and by formulating the equation in a conservative form [36]. The approach presented in this paper is based on theconservative phase field approach by Chiu and Lin [17], which incorporates an anti-diffusive term to Eq.1, yielding ∂φ∂t + u · ∇ φ = γ (cid:20) ∇ φ − φ (1 − φ )(1 − φ )2 (cid:15) − |∇ φ |∇ · n (cid:21) . (2)where n = ∇ φ |∇ φ | is the interface’s unit normal vector. Using a kernel to define the phase field in terms of a distance function, ψ such that φ = (cid:2) tanh (cid:0) ψ (cid:15) (cid:1)(cid:3) , the first and second order derivatives can be expressed algebraically as follows |∇ φ | = ∂φ∂ψ = φ (1 − φ ) (cid:15) (3) |∇ φ | = ∂ φ∂ψ = φ (1 − φ )(1 − φ ) (cid:15) (4)by collecting the terms and assuming a divergence-free velocity field, the following equation is derived ∂φ∂t + ∇ · ( φ u ) = γ ( (cid:15) ∇ · ( ∇ φ ) − ∇ · [ φ (1 − φ ) n ]) (5)where γ = γ(cid:15) , γ = | u | max . In this work, however, it was found that a local choice of the multiplier yields better results,especially in areas of the flow with very small velocity magnitudes, rewriting Eq. 5 as ∂φ∂t + ∇ · ( φ u ) = max( γ, | u | ) ( (cid:15) ∇ · ( ∇ φ ) − γ ∇ · [ φ (1 − φ ) n ]) . (6)Numerical experiments showed that the optimal parameters are γ = 0 .
09 and (cid:15) = 0 . h/p where h is the average grid size and p is the order of the nodal basis used in the FR discretization, which will be discussed in the following sections.A major issue with Eq. 5 is the ill-conditioning of the normal vector, especially away from the interface, which causeserroneous estimate of interface location and thickness. This issue is common with the conservative level-set method [20],due to the similarity of this equation with the re-initialization equation used to maintain the hyperbolic tangent profile ofthe interface. In the conservative level-set literature, attempts to improve the accuracy of normal computations were madeby recycling the signed distance function using inverse hyperbolic tangent [37], restoring the signed distance function by afast marching method (FMM) [38], and algebraic stabilization terms [20, 39]. Those approaches, however, do not completelyeliminate spurious and oscillatory interface behavior.Inspired by the notion of non-linear pre-conditioning [40], Chiu[41] replaced the right hand side of Eq. 5 with an equivalentexpression in terms of the distance function in addition to a space-time Lagrange multiplier. To provide ψ , this approachrequired an additional phase field equation.In this paper, we propose an alternative preconditioning method based on traditional level-set re-initialization which minimizesinterface displacement and does away with the need to compute expensive Lagrangian multipliers. The normal vector in Eq.6 is instead computed using ψ , such that n = ∇ ψ |∇ ψ | (7)and ψ is simply advected according to ∂ψ∂t + ∇ · ( ψ u ) = 0 . (8)In the beginning of the simulation, and every N time-steps, ψ is reset using ψ = φ −
12 (9)to form the initial condition of the following stabilized preconditioning equation ∂ψ∂τ + sgn ( ψ )(1 − |∇ ψ | ) = ∇ · ( ν s ∇ ψ ) . (10) here τ is pseudo-time and is irrelevant to the physical problem, sgn ( ψ ) = tanh( ψ |∇ ψ | ) is a smeared sign function and ν s isthe stabilization viscosity. Eq. 10 is similar in form to the level-set re-initialization equation, however, it is not iterated untilconvergence to the signed-distance field. Instead we found it sufficient to only carry out the preconditioning procedure every500 − We start by rewriting the phase field time-evolution (Eq. 6) and preconditioning (Eq. 10) equations in the following form ∂ u ∂t + ∇ · f ( u , q ) = 0 (11) q − ∇ u = 0 (12)where u = [ φ, ψ ] T . Borrowing an analogy to the Navier-Stokes discretization, the fluxes of the phase field and preconditioningvariables are split into ”inviscid” and ”viscous” parts f = f ( inv ) − f ( vis ) , such that f ( inv ) = (cid:20) φv x φv y φv z ψv x ψv y ψv z (cid:21) , (13) f ( vis ) = max( γ, | v | ) (cid:34) (cid:15) ∂φ∂x − γ (cid:104) φ (1 − φ ) |∇ ψ | ∂ψ∂x (cid:105) (cid:15) ∂φ∂y − γ (cid:104) φ (1 − φ ) |∇ ψ | ∂ψ∂y (cid:105) (cid:15) ∂φ∂z − γ (cid:104) φ (1 − φ ) |∇ ψ | ∂ψ∂z (cid:105) (cid:35) . (14)Let the physical domain be divided into | Ω | elements, where Ω is the set of non-overlapping elements, such thatΩ = | Ω |− (cid:91) n =0 Ω n and | Ω |− (cid:92) n =0 Ω n = ∅ (15)where for the sake of brevity we assume that Ω n can be of any geometric shape in the dimension of the problem, such asquadrilateral or triangle in 2D, and prism, hexahedron, etc. in 3D.The partial differential equations are satisfied inside each element, for instance ∂u n ∂t + ∇ · f ( u n , q n ) = 0 , (16) q n − ∇ u n = 0 (17)where from this point on, we use the non-bold typface u to refer to either φ or ψ . In order to simplify implementation andimprove computational efficiency, flux reconstruction steps are carried out in transformed element space. This is achieved bymeans of mapping functions that transform different element types into their respective standard elements. These standardelements are defined in transformed coordinates, (cid:101) x , such that (cid:101) x = M − ( x ) (18) x = M ( (cid:101) x ) (19)The Jacobians related to this mapping are J n = ∂ M ni ∂ (cid:101) x j , J n = det( J n ) (20) J − n = ∂ M − ni ∂ (cid:101) x j , J − n = det( J − n ) = 1 / J n . (21)Using the transformation above, the flux and auxiliary variable in transformed space are (cid:101) f n ( (cid:101) x , t ) = J n ( (cid:101) x ) J − n ( M n ( (cid:101) x )) f n ( M n ( (cid:101) x ) , t ) (22) (cid:101) q n ( (cid:101) x , t ) = J Tn ( (cid:101) x ) ∇ u n ( M n ( (cid:101) x ) , t ) . (23) sing (cid:101) ∇ = ∂∂ (cid:101) x i , the physical solution’s time-derivative in Eq. 16 can be expressed in terms of the transformed divergenceof the flux in transformed, standard element space ∂u n ∂t + J − n (cid:101) ∇ · (cid:101) f n ( (cid:101) x , t ) . (24)In the flux reconstruction method, a set of solution points, (cid:101) x ( u ) ni , where 1 < i < N ( u ) n , are placed within standard elementsusing an appropriate distribution. In this work, points are distributed according to the Gauss-Legendre points in quadrilateralelements and Wiallms-Shunn points in triangles. Next, a nodal basis set l ni ( (cid:101) x ) is constructed using Ψ( (cid:101) x ) nj , which is a basisset that spans a polynomial space of order p such that l ni ( (cid:101) x ) = V − nij Ψ nj ( (cid:101) x ) (25)where V nij = Ψ ni ( (cid:101) x ( u ) nj ) are the elements of the Vandermonde matrix, and are required to satisfy the property l ni ( (cid:101) x nj ) = δ ij .In addition to solution points, a set of flux points , (cid:101) x ( f ) ni where 1 < i < N ( f ) n , are defined on element boundaries. Fluxpoints are distributed in transformed space in such a way that they share the same physical coordinates with the correspondingpoints of neighboring elements.The steps of solving equations 22 and 23 are as follows: Step 1:
In order to compute gradients (Eq. 23), the continuous solution across element boundaries is first reconstructed.The discontinuous solution in the n th element at the i th flux point u ( f ) ni is approximated by interpolating from the discontinuoussolution at solution points u ( u ) nj using the polynomial nodal basis u ( f ) ni = N ( f ) n (cid:88) j =1 u ( u ) nj l nj ( (cid:101) x ( f ) ni ) (26) Step 2:
By adopting the Local Discontinuous Galerkin approach, the common solution at flux points is then computed C ( u ( f ) ni ) = C ( u ( f ) n (cid:48) i (cid:48) ) = ( 12 − β ) u ( f ) ni + ( 12 + β ) u ( f ) n (cid:48) i (cid:48) (27)where n (cid:48) i (cid:48) denote the coinciding flux point of a neighboring element, and the parameter β = ± Step 3:
The gradient of the continuous solution at solution points is computed using (cid:101) q ( u ) ni = N ( f ) n (cid:88) j =1 (cid:104)(cid:98)(cid:101) n nj · (cid:101) ∇ · g nj ( ˜x ( u ) ni ) (cid:16) C ( u ( f ) nj ) − u ( f ) nj (cid:17)(cid:105) + N ( u ) n (cid:88) j =1 (cid:104) u ( u ) nj (cid:101) ∇ l ( u ) nj ( ˜x ( u ) ni ) (cid:105) (28)where (cid:101) n is the normalized and transformed outward-pointing normal vector at the flux points, and g nj ( x ( u ) ) is the vectorcorrection function associated with each flux point, which is used to ensure gradient continuity across element interfaces andis required to satisfy (cid:98)(cid:101) n ni · g nj ( x ( f ) ni ) = δ ij (29)The transformed solution’s gradient is mapped back to physical space using q ( u ) ni = J − T (cid:101) q ( u ) ni , and then evaluated at fluxpoints, using the same procedure in 26 q ( f ) ni = N ( f ) n (cid:88) j =0 q ( u ) nj l nj ( (cid:101) x ( f ) ni ) (30) Step 4:
Using the transformed gradient found in the previous steps, the flux is evaluated at solution points in transformedspace according to (cid:101) f ( u ) ni = J ni J − ni f ( u ( u ) ni , q ( u ) ni ) . (31)Using this, the normal transformed flux at flux points is computed in a similar manner to Eqs. 26 and 30 (cid:101) f f ⊥ ni = N ( f ) n (cid:88) j =0 (cid:98)(cid:101) n · (cid:101) f ( u ) ni l nj ( (cid:101) x ( f ) ni ) (32) tep 5: In order to find the continuous flux, common normal flux at flux points is computed using F ( f f ⊥ ni ) = − F ( f f ⊥ n (cid:48) i (cid:48) ) = F ( inv ) ( f ( f )( inv ) ni ) − F ( vis ) ( f ( f )( vis ) ni ) (33)Note that unlike the discontinuous normal flux found in Eq. 32, the common normal flux is found using quantitiesextrapolated from solution points from both sides of the interface, i.e. u ( f ) in , u ( f ) i (cid:48) n (cid:48) , q ( f ) in and q ( f ) i (cid:48) n (cid:48) . The inviscid part is simplythe average of inviscid normal fluxes F ( inv ) ( f ( f )( inv ) ni ) = (cid:98) n ni · (cid:34) f ( f )( inv ) ni + f ( f )( inv ) n (cid:48) i (cid:48) (cid:35) (34)and the viscous part is found using the LDG approach F ( vis ) ( f ( f )( vis ) ni ) = (cid:98) n · (cid:20)(cid:18)
12 + β (cid:19) f ( f ) ni + (cid:18) − β (cid:19) f ( f ) n (cid:48) i (cid:48) (cid:21) + τ ( u ( f ) ni − u ( f ) n (cid:48) i (cid:48) ) . (35)Common normal fluxes are then scaled and transformed to standard element space to prepare for the final step. Step 6:
The divergence of the continuous flux is obtained via a procedure that is analogous to Eq. 28 (cid:101) ∇ · ˜f ( u ) ni = N ( f ) n (cid:88) j =0 (cid:104) (cid:101) ∇ · ˜g nj ( x ( u ) ni ) (cid:16) F ( ˜ f f ⊥ nj ) − ˜ f ( f ⊥ ) nj (cid:17)(cid:105) + N ( u ) n (cid:88) j =0 (cid:104) ˜f ( u ) nj · (cid:101) ∇ l ( u ) nj ( ˜x ( u ) ni ) (cid:105) (36)the divergence of the flux at each solution point is then marched in time using a 4 stage Runge-Kutta scheme ∂u ( u ) ni ∂t = −J − ni (cid:101) ∇ · (cid:101) f ( u ) ni (37) The system of equations used for preconditioning reads ∂ψ∂τ + ∇ · g ( ψ, q ) − S ( ψ, q ) = 0 (38) q − ∇ ψ = 0 (39)where g = g ( vis ) = − ν s ∇ ψ , and S ( u , q ) = sgn ( ψ )(1 − | q | ). The steps taken to reconstruct the continuous ψ field andcomputing its gradient are similar to those for the phase field equation (steps 1 -3). However, step 4 is preceded by theintermediate step of computing the element-specific artificial viscosity amount needed to stabilize the equation by suppressinghigh-frequency oscillations. These oscillations arise due to the discontinuity of q and source term S ( ψ, q ) in Eq. 38. The jumpin ψ across element boundaries is exacerbated with pseudo-time marching, especially in regions with very small gradients orwhenever kinks arise. To avoid excessive diffusion in regions where it is not needed, we use localized artificial viscosity forstabilization, which is used in PyFR to treat discontinuities arising due to shocks in high Mach number compressible flows. Abrief summary of the method by for finding ν s is given here, but more details can be found in [42]. The basic idea behind theshock sensor is to quantify the decay rate of orthogonal basis expansion coefficients. In smooth fields, these coefficients decayrapidly, but when the solution is under-resolved, posing as a discontinuity, the coefficients decay in relation to the strength ofthe discontinuity. The resolution indicator is defined as follows s e = log (cid:18) (cid:104) u − ˆ u, u − ˆ u (cid:105)(cid:104) u, u (cid:105) (cid:19) (40)where (cid:104)· , ·(cid:105) is the inner product, u is the solution described in modal basis of order N , and ˆ u is the truncated solution,projected onto the polynomial space of order N −
1, as follows u = N (cid:88) i =1 χ i Ψ i (41)ˆ u = N − (cid:88) i =1 χ i Ψ i (42) χ i being the modal expansion coefficients, found using the definition of the Vandermonde matrix, u = V χ . smoothed, element-wise stabilization viscosity is then found according to ν s = if s e < s − κ ν max (cid:16) π ( s e − s )2 κ (cid:17) if s − κ ≤ s e ≤ s + κν max if s e > s + κ (43)where we set s = − κ = 5 × and ν max = 0 . h/p After completing the sixth stage of the flux reconstruction, the semi-discrete form of the pseudo-time evolution of thepreconditioning equation is obtained ∂ψ ( u ) ni ∂τ = −J − ni (cid:101) ∇ · ˜ g ni + S ni . (44)We found it sufficient to only carry out the above preconditioning procedure every 1000 physical time-steps and iteratingEq. 10 until the L norm satisfies (cid:118)(cid:117)(cid:117)(cid:116) (cid:80) | Ω | n =1 (cid:80) N ( u ) n i =1 (cid:16) R ( u ) ni (cid:17) N ( u ) ≤
12 (45)where the residual is R ( u ) ni = ψ n − ψ n − ∆ τ (46)and n is the pseudo-time iteration number, which typically requires 5 −
20 iterations using the four stage Runge Kuttascheme.
In order to evaluate the common solution C ( u ( f ) ni ) and common normal flux F ( f f ⊥ ni ) at the nodes of boundary faces, ghoststates are used to define the solution B ( u ( f ) ni ) and its gradient B ( q ( f ) ni ) on virtual nodes on the other side of the interface. Forall simulations shown in this paper, the following ghost states were used for the phase field system B ( u ( f ) ni ) = u ( f ) ni (47) B ( q ( f ) ni ) = q ( f ) ni (48) In this section, the proposed method is validated by comparing to popular 2D and 3D interface capturing benchmarks andthe results are compared to previous studies. In order to provide a quantitative error measure that takes into account theconservativeness of our approach and penalizes the presence of artifacts such as smearing and fragmentation, we use thefollowing error measures. The L error, is computed as follows E L = (cid:90) | φ i − φ e | d Ω (49)where φ i and φ e are the numerical and exact solutions respectively. Another measure that is widely used when quantifyinginterface error is the relative error E r = (cid:82) | φ i − φ e | d Ω (cid:82) φ e d Ω . (50)To quantify mass conservation, we use the following measure E m = (cid:82) φ i d Ω − (cid:82) φ d Ω (cid:82) φ e d Ω (51)where φ is the phase field at the initial condition. Using the above-mentioned error measures also allows comparison toa large number of previously published interface capturing works.The phase field at solution points is initialized according to φ i = 12 (1 + tanh( βψ di )) (52) here ψ d is the signed distance function of the initial condition and β is an interface sharpness parameter, which is set to β = 3 p/ h in order to produce a sharp yet smooth and non-oscillatory initial condition. Zalesak’s test problem [43] is used in this section to assess the accuracy and long-term time integration behavior of the proposedmethod. At the initial condition, a slotted disk with a radius r = 0 .
15 centered at (0 . , .
75) in a unit square domain. Thedisk is initialized using Eq. 52 with the following distance function ψ di = r − (cid:112) ( x − x ) + ( y − y ) . (53)The slot is defined by ( | x − . | ≤ . ∧ ( y ≤ . . − y, x − .
5) until it returns to its original position. This test was carried out on 100 ×
100 and 200 ×
200 uniformquadrilateral meshes.Fig. 1 compares the initial and final shapes of Zalesak’s disk on the 100 mesh using p = 2 and p = 4. A qualitativeexamination of the φ = 0 . p = 2, but some deviation from the initial conditioncan be seen in the corners of the slot. Increasing the basis order to p = 4 significantly improves the ability to resolve suchchallenging features. Results using finer meshes or higher basis orders are not shown here as the initial and final interfacecontours become visually indistinguishable.The effect of long-time integration on challenging interface features is investigated by examining the distortion in thecorners of the slot after multiple rotations in Fig. 2. This demonstrates that even with long-term time integration, ourapproach is able to preserve such challenges features and maintain a constant interface thickness.Numerical errors using different mesh resolutions and orders are provided in Table 1 and compared to previous studies.Despite its simplicity, the FR-PCPF approach produces more accurate results when compared to more complicated approachesutilizing high order polynomial interface reconstruction.Figures 3 and 4 show the phase field, φ and auxiliary preconditioning variable ψ plotted along the center-lines of the slotteddisk at x = 0 . y = 0 .
75. The preconditioning procedure produces a sufficiently smooth ψ field which allows for accuratenormal vector computation with minimal offset between the interfaces of the phase field and auxiliary variable at φ = 0 . ψ = 0 respectively.We also carry out the slotted disk problem proposed by Rudman [44]. In this problem, a slotted disk with r = 0 . . , .
75) with a slot defined by ( | x − . | ≤ . ∧ ( y ≤ .
65) is placed in a [0 , domain with a 200 ×
200 uniformquadrilateral mesh. The velocity field in this case is (2 − y, x − p = 2 and p = 4, which again demonstrate the accuracy of the FR-PCPF interface capturingapproach and its ability to resolve challenging interface features. This is further confirmed by quantitatively examining theerror in comparison to published literature as shown in Table 2. Our results using p = 2 outperform geometric VOF methodswith a significant margin, while producing comparable results with the coupled THINC/level set approach [19] method withsecond order polynomial interface reconstruction. Using p ≥ Method MTHINC [45] 1 . × − . × − THINC/LS(P2) [19] 1 . × − . × − THINC/LS(P4)[19] 9 . × − . × − FR-PCPF( p = 2) 1 . × − . × − FR-PCPF( p = 3) 3 . × − . × − FR-PCPF( p = 4) 2 . × − . × − FR-PCPF( p = 5) 1 . × − . × − Table 1: E r error for the Zalesak problem on a quadrilateral mesh after one rotation8 0 . . . . . . . . Resolution : 100 , p = 2 . . . . . . . . Resolution : 100 , p = 4 Figure 1: Initial (black) versus final (red) φ = 0 . .
47 0 .
48 0 .
49 0 . .
51 0 .
52 0 . . . . . .
88 0 .
51 0 .
52 0 .
53 0 .
54 0 .
55 0 . . . . . . . .
05 (dashed), 0 . .
95 (dash-dotted) φ contours of the initial condition(black) of Zalesak’s disk on the 100 mesh and using P = 4 and after 1(blue), 2(red), 3(green), 4(magenta) rotations.9 .0 0.2 0.4 0.6 0.8 1.00.00.20.40.60.81.0 x φ -0.5-0.3-0.10.10.30.5 ψ Figure 3: φ and ψ plotted along the x-axis of Zalesak’s disk after 1 rotation φ -0.5-0.3-0.10.10.30.5 ψ Figure 4: φ and ψ plotted along the y-axis of Zalesak’s disk after 1 rotation10 0 . . . . . . . . Resolution : 200 , p = 2 . . . . . . . . Resolution : 200 , p = 4 Figure 5: Initial (black) versus final (red) φ = 0 . Method
ErrorEMFPA/Youngs [46] 1 . × − EMFPA/Puckett [46] 9 . × − EMFPA/SIR [46] 8 . × − Geometric predictor-corrector [47] 9 . × − THINC/QQ [48] 1 . × − Quadratic fit + Lagrangian[49] 4 . × − THINC/LS(P2) [19] 4 . × − THINC/LS(P4) [19] 3 . × − FR-PCPF( p = 2) 4 . × − FR-PCPF( p = 3) 1 . × − FR-PCPF( p = 4) 1 . × − FR-PCPF( p = 5) 8 . × − Table 2: E r for the Rudman problem on a 200 ×
200 quadrilateral mesh after one rotation
In this benchmark case, which was first proposed in [50], a disk with a radius r = 0 .
15 centered at (0 . , .
75) is initially placedwithin a domain of size [0 , . The disk is then advected by the following velocity field u ( x, y, t ) = − sin ( πx ) sin(2 πy ) cos( πtT ) v ( x, y, t ) = sin(2 πx ) sin ( πy ) cos( πtT ) . (54)The spatially and temporally varying velocity field described above significantly deforms the disk by stretching it into a thinspiral filament until t = T /
2. In this second half of the test the flow is reversed, such that if a perfect scheme were to beused the disk would recover its initial shape at t = T . In all simulations shown in this subsection the period was set to = 8. Low-order methods are typically unable to resolve sub-grid interface features, thereby producing fragmented dropletswhenever the width of the filament drops below that of the grid. However, since in the present approach all field variables areprojected to a high-order polynomial space, sub-grid interface features can be resolved and such artifacts are not produced,provided that a sufficiently high order basis is used.Fig. 7, illustrates the significant improvement in accuracy using the preconditioned phase field method, when compared tothe original method. The large variations in the thickness of the interface that result from using the original formulation areeliminated due to preconditioning. Furthermore, the proposed method also eliminates high-frequency oscillation, gaps, andfragmentation in the phase field.Three meshes shown in Fig. 6 were used for this test to demonstrate the robustness of this approach against mesh type.The three types are: • Uniform Quadrilateral Mesh (Fig. 6a): it is a uniform Cartesian mesh, which has been thoroughly investigated inliterature and is used to provide a baseline for comparison with other mesh types. • Nearly Equilateral Triangular Mesh (Fig. 6b): it is an unstructured triangular mesh generated using the advancing-frontalgorithm. • Uniform Split Triangular Mesh (Fig. 6c): it is an unstructured triangular mesh produced by diagonally splitting aCartesian mesh.Errors for different meshes and basis orders are listed in Table 3, quantitatively demonstrating the improvement in accuracyresulting from our proposed preconditioning method. The improvement can range from 80% to over 300% in some cases. Thetable shows that such improvement does not seem to be consistent with increasing order or base mesh resolution and couldrequire further investigation.Fig. 8 illustrates the improvement of the accuracy with increased mesh resolution and basis orders up to 64 , p = 4 as nofurther improvement can be discerned visually with finer meshes or higher orders.Table 4 compares the E L error obtained from the proposed scheme to some of the best and most recent interface capturingstudies. Using p = 2 is adequate to achieve comparable or superior accuracy when compared to other methods. Consistentimprovement in accuracy can be seen when the polynomial order is increased.In Fig. 9 the ability of the FR-PCPF method to resolve sub-grid interface features is demonstrated by providing an enlargedview of the tip of the vortex at t = T /
2. Multiple interfaces within an element are resolved without introducing artifacts intothe solution. (a) (b) (c)
Figure 6: Meshes used in Rider-Koth cases
The time history of global mass change ( E m ) in the domain is provided in figs. 11 and 11. The figures clearly demonstratesimprovement in global mass conservation with increased effective resolution ( h/p ). In order to verify the ability of the proposed method to capture distorted interfaces in 3D, the deformation of a sphere in a3D velocity field is simulated [53]. A sphere with radius r = 0 .
15 in a unit domain centered at (0 . , . , .
35) is deformed
12 0 . . . . . . . . . . . . . . . . φ = 0 . , . , .
95 contour lines at t = T for the Rider-Kothevortex problem Preconditioned CPF Unpreconditioned CPF Mesh p = 3 p = 4 p = 5 p = 3 p = 4 p = 5 b − . × − . × − . × − . × − . × − . × − b − . × − . × − . × − . × − . × − . × − b − . × − . × − . × − . × − . × − . × − c − . × − . × − . × − . × − . × − . × − c − . × − . × − . × − . × − . × − . × − c − . × − . × − . × − . × − . × − . × − Table 3: E L error for the Rider-Koth vortex on different meshes, resolutions and polynomial orders using thepreconditioned and unpreconditioned conservative phase field equation Method Rider-Kothe/Puckett [50] 4 . × − . × − . × − EMFPA/Puckett [51] 3 . × − . × − . × − THINC/QQ [48] 6 . × − . × − . × − Markers-VOF[52] 7 . × − . × − . × − ISLSVOF [18] 4 . × − . × − . × − THINC/LS (P2) [19] 1 . × − . × − . × − THINC/LS (P4) [19] 2 . × − . × − . × − FR-PCPF( p = 2) 7 . × − . × − . × − FR-PCPF( p = 3) 7 . × − . × − . × − FR-PCPF( p = 4) 5 . × − . × − . × − FR-PCPF( p = 5) 3 . × − . × − . × − Table 4: E L error for the Rider-Koth vortex test on quadrilateral mesh13 0 . . . . . . . . Resolution : 32 , p = 2 . . . . . . . . Resolution : 32 , p = 4 . . . . . . . . Resolution : 64 , p = 2 . . . . . . . . Resolution : 64 , p = 4 Figure 8: Final (red) and initial (black) φ = 0 . φ = 0 . , p = 5 case with an enalrged view of the tail0 0 . . . . . . . . . − · − t/T E m Figure 10: Mass error on a quad 64 mesh for p = 3(black), p = 4(blue), p = 5(red)0 0 . . . . . . . . . − · − t/T E m Figure 11: Mass error on 32 (black), 64 (blue), 128 (red) quad meshes using p = 415 ith the following velocity field u ( x, y, t ) = 2 sin ( πx ) sin(2 πy ) sin(2 πz ) cos( πtT ) v ( x, y, t ) = − sin(2 πx ) sin ( πy ) sin(2 πz ) cos( πtT ) w ( x, y, t ) = − sin(2 πx ) sin(2 πy ) sin ( πz ) cos( πtT ) . (55)For all tests in this section, the period T = 3. Similarly to the Rider-Koth vortex, the sphere is deformed until T / t = T . This test was done on hexahedral mesheswith 32 , 64 and 128 elements using p − p
5. The large number of solution points needed to carry out these tests (up to0 .
45 billion points for the 128 mesh with p = 5 case) necessitated using as many as 40 Nvidia P100 GPUs to accommodatethe large computational resources needed.The evolution of the sphere for the 128 mesh case and p = 4 is depicted in Fig. 12. The figure shows some minor artifactsat the interface while remaining quantitatively accurate as highlighted by the values of E L error comparison outlined inTable 5. The table further shows that the FR-PCPF method is able to produce more accurate results even in comparisonto complex methods that rely on geometric reconstruction or employ high-order schemes that are also able resolve sub-cellinterface structures. Method LVIRA [54] 6 . × − . × − . × − THINC/QQ [48] 7 . × − . × − . × − DS-CLSMOF [55] 4 . × − . × − − ISLSVOF [18] 8 . × − . × − . × − THINC/LS (P2) [19] 6 . × − . × − . × − THINC/LS (P4) [19] 5 . × − . × − . × − FR-PCPF( p = 2) 6 . × − . × − . × − FR-PCPF( p = 3) 3 . × − . × − . × − FR-PCPF( p = 4) 2 . × − . × − . × − FR-PCPF( p = 5) 1 . × − . × − . × − Table 5: E L error for the 3D vortex test on hexahedral meshes To pave the way to high-order multi-phase simulations using the Flux Reconstruction method, we develop a novel interfacecapturing approach. The merits and novelties of the proposed method can be summarized as follows: • To enable capturing sharp oscillation-free interfaces using a high-order approach, we develop a localized conservativephase field method. • The accuracy of the phase field method is significantly improved by introducing a level-set based preconditioning stepthat increases the accuracy of interface normals, and prevents the appearance of smearing and fragmentation artifacts. • The favorable numerical properties of the FR-PCPF approach, such as the high-order polynomial representation of allfield variables, makes it possible to capture sub-grid interface features with a high level of accuracy.Numerical validation test cases have shown that using second order polynomial basis is sufficient to produce results thatare comparable to some of the most accurate existing methods. Increasing the polynomial order to higher than 3 consistentlyimproves numerical accuracy. The FR-PCPF method is shown to be able to maintain good mass conservation and resolvechallenging interface features even with long time integration. Increasing effective resolution ( h/p ) was shown to produceconsistent improvements in global mass conservation, accuracy and the ability to resolve complicated interface features. =0.0t=0.5 t=1.0 t=1.5 t=2.0 t=2.5t=3 Figure 12: Time evolution of the φ = 0 . eferences [1] F. D. Witherden, A. M. Farrington, and P. E. Vincent, “PyFR : An open source framework for solving advection– diffusion type problems on streaming architectures using the flux reconstruction approach ,” Computer PhysicsCommunications , vol. 185, no. 11, pp. 3028–3040, 2014. [Online]. Available: http://dx.doi.org/10.1016/j.cpc.2014.07.011[2] F. Harlow and J. Welch, “Volume tracking methods for interfacial flow calculations,”
Physics of fluids , vol. 8, pp. 21–82,1965.[3] C. Hirt, A. Amsden, and J. Cook, “An arbitrary lagrangian–eulerian computing method for all flow speeds,”
Journal ofcomputational physics , vol. 135, no. 2, pp. 203–216, 1997.[4] J. Baiges, R. Codina, A. Pont, and E. Castillo, “An adaptive fixed-mesh ale method for free surface flows,”
ComputerMethods in Applied Mechanics and Engineering , vol. 313, pp. 159–188, 2017.[5] P. Nithiarasu, “An arbitrary lagrangian eulerian (ale) formulation for free surface flows using the characteristic-basedsplit (cbs) scheme,”
International Journal for Numerical Methods in Fluids , vol. 48, no. 12, pp. 1415–1428, 2005.[6] E. Marchandise, J. F. Remacle, and N. Chevaugeon, “A quadrature-free discontinuous Galerkin method for the level setequation,”
Journal of Computational Physics , vol. 212, no. 1, pp. 338–357, 2006.[7] R. L¨ohner, C. Yang, and E. O˜nate, “On the simulation of flows with violent free surface motion,”
Computer Methods inApplied Mechanics and Engineering , vol. 195, no. 41-43, pp. 5597–5620, 2006.[8] K. Kleefsman, G. Fekken, A. Veldman, B. Iwanowski, and B. Buchner, “A volume-of-fluid based simulation method forwave impact problems,”
Journal of computational physics , vol. 206, no. 1, pp. 363–393, 2005.[9] M. M. Kamra, N. Mohd, C. Liu, M. Sueyoshi, and C. Hu, “Numerical and experimental investigation of three-dimensionality in the dam-break flow against a vertical wall,”
Journal of Hydrodynamics , vol. 30, no. 4, pp. 682–693,2018.[10] R. P. Fedkiw, T. Aslam, and S. Xu, “The ghost fluid method for deflagration and detonation discontinuities,”
Journal ofComputational Physics , vol. 154, no. 2, pp. 393–427, 1999.[11] Y. Pei, S. Som, E. Pomraning, P. K. Senecal, S. A. Skeen, J. Manin, and L. M. Pickett, “Large eddy simulation of areacting spray flame with multiple realizations under compression ignition engine conditions,”
Combustion and Flame ,vol. 162, no. 12, pp. 4442–4455, 2015.[12] K. Luo, C. Shao, M. Chai, and J. Fan, “Level set method for atomization and evaporation simulations,”
Progress inEnergy and Combustion Science , vol. 73, pp. 65–94, 2019. [Online]. Available: https://doi.org/10.1016/j.pecs.2019.03.001[13] W. Yang, M. Jia, K. Sun, and T. Wang, “Influence of density ratio on the secondary atomization of liquid droplets underhighly unstable conditions,”
Fuel , vol. 174, pp. 25–35, 2016.[14] F. Xiao, Y. Honma, and T. Kono, “A simple algebraic interface capturing scheme using hyperbolic tangent function,”
International Journal for Numerical Methods in Fluids , vol. 48, no. 9, pp. 1023–1040, 2005.[15] E. Olsson and G. Kreiss, “A conservative level set method for two phase flow,”
Journal of Computational Physics , vol.210, no. 1, pp. 225–246, 2005.[16] Y. Sun and C. Beckermann, “Sharp interface tracking using the phase-field equation,”
Journal of Computational Physics ,vol. 220, no. 2, pp. 626–653, 2007.[17] P. H. Chiu and Y. T. Lin, “A conservative phase field method for solving incompressible two-phase flows,”
Journal ofComputational Physics , vol. 230, no. 1, pp. 185–204, 2011. [Online]. Available: http://dx.doi.org/10.1016/j.jcp.2010.09.021[18] K. G. Lyras, B. Hanson, M. Fairweather, and P. J. Heggs, “A coupled level set and volume of fluid method with are-initialisation step suitable for unstructured meshes,”
Journal of Computational Physics , vol. 407, p. 109224, 2020.[Online]. Available: https://doi.org/10.1016/j.jcp.2019.109224[19] L. Qian, Y. Wei, and F. Xiao, “Coupled THINC and level set method: A conservative interface capturing scheme withhigh-order surface representations,”
Journal of Computational Physics , vol. 373, pp. 284–303, 2018. [Online]. Available:https://doi.org/10.1016/j.jcp.2018.06.074[20] N. Shervani-Tabar and O. V. Vasilyev, “Stabilized conservative level set method,”
Journal of Computational Physics ,vol. 375, pp. 1033–1044, 2018. [Online]. Available: https://doi.org/10.1016/j.jcp.2018.09.020[21] A. Karakus, T. Warburton, M. H. Aksel, and C. Sert, “A GPU accelerated level set reinitialization for an adaptivediscontinuous Galerkin method,”
Computers and Mathematics with Applications , vol. 72, no. 3, pp. 755–767, 2016.[Online]. Available: http://dx.doi.org/10.1016/j.camwa.2016.05.025[22] P. E. Vincent, P. Castonguay, and A. Jameson, “A new class of high-order energy stable flux reconstruction schemes,”
Journal of Scientific Computing , vol. 47, no. 1, pp. 50–72, 2011.
23] L. Wang and M. Yu, “A comparative study of implicit jacobian-free rosenbrock-wanner, esdirk and bdf methods forunsteady flow simulation with high-order flux reconstruction formulations,” arXiv preprint arXiv:1904.04825 , 2019.[24] S. Matsushita and T. Aoki, “A weakly compressible scheme with a diffuse-interface method for low Machnumber two-phase flows,”
Journal of Computational Physics , vol. 376, pp. 838–862, 2019. [Online]. Available:https://doi.org/10.1016/j.jcp.2018.10.019[25] Z. Jibben and M. Herrmann, “An arbitrary-order Runge–Kutta discontinuous Galerkin approach to reinitialization forbanded conservative level sets,”
Journal of Computational Physics , vol. 349, pp. 453–473, 2017. [Online]. Available:http://dx.doi.org/10.1016/j.jcp.2017.08.035[26] J. Zhang and P. Yue, “A high-order and interface-preserving discontinuous Galerkin method for level-set reinitialization,”
Journal of Computational Physics , vol. 378, pp. 634–664, 2019. [Online]. Available:https://doi.org/10.1016/j.jcp.2018.11.029[27] A. Karakus, T. Warburton, M. H. Aksel, and C. Sert, “An adaptive fully discontinuous Galerkin level set method forincompressible multiphase flows,”
International Journal of Numerical Methods for Heat and Fluid Flow , vol. 28, no. 6,pp. 1256–1278, 2018.[28] J. Grooss and J. S. Hesthaven, “A level set discontinuous Galerkin method for free surface flows,”
Computer Methods inApplied Mechanics and Engineering , vol. 195, no. 25-28, pp. 3406–3429, 2006.[29] D. Gottlieb and C.-W. Shu, “On the gibbs phenomenon and its resolution,”
SIAM review , vol. 39, no. 4, pp. 644–668,1997.[30] H. T. Huynh, “A flux reconstruction approach to high-order schemes including discontinuous galerkin methods,” in , 2007, p. 4079.[31] N. A. Loppi, F. D. Witherden, A. Jameson, and P. E. Vincent, “A high-order cross-platform incompressible Navier –Stokes solver via artificial compressibility with application to a turbulent jet,”
Computer Physics Communications , vol.233, pp. 193–205, 2018. [Online]. Available: https://doi.org/10.1016/j.cpc.2018.06.016[32] S. M. Allen and J. W. Cahn, “A microscopic theory for antiphase boundary motion and its application to antiphasedomain coarsening,”
Acta metallurgica , vol. 27, no. 6, pp. 1085–1095, 1979.[33] J. W. Cahn and J. E. Hilliard, “Free energy of a nonuniform system. i. interfacial free energy,”
The Journal of chemicalphysics , vol. 28, no. 2, pp. 258–267, 1958.[34] J. Rubinstein and P. Sternberg, “Nonlocal reaction—diffusion equations and nucleation,”
IMA Journal of Applied Math-ematics , vol. 48, no. 3, pp. 249–264, 1992.[35] M. Brassel and E. Bretin, “A modified phase field approximation for mean curvature flow with conservation of thevolume,”
Mathematical Methods in the Applied Sciences , vol. 10, no. 34, pp. 1157–1180, 2011.[36] D. Jeong and J. Kim, “Conservative allen–cahn–navier–stokes system for incompressible two-phase fluid flows,”
Computers& Fluids , vol. 156, pp. 239–246, 2017.[37] L. Zhao, X. Bai, T. Li, and J. Williams, “Improved conservative level set method,”
International Journal for NumericalMethods in Fluids , vol. 75, no. 8, pp. 575–590, 2014.[38] O. Desjardins, V. Moureau, and H. Pitsch, “An accurate conservative level set/ghost fluid method for simulating turbulentatomization,”
Journal of Computational Physics , vol. 227, no. 18, pp. 8395–8416, 2008.[39] R. K. Shukla, C. Pantano, and J. B. Freund, “An interface capturing method for the simulation of multi-phase compressibleflows,”
Journal of Computational Physics , vol. 229, no. 19, pp. 7411–7439, 2010.[40] K. Glasner, “Nonlinear preconditioning for diffuse interfaces,”
Journal of Computational Physics , vol. 174, no. 2, pp.695–711, 2001.[41] P. H. Chiu, “A coupled phase field framework for solving incompressible two-phase flows,”
Journal of ComputationalPhysics , vol. 392, pp. 115–140, 2019. [Online]. Available: https://doi.org/10.1016/j.jcp.2019.04.069[42] P. O. Persson, “Shock capturing for high-order discontinuous Galerkin simulation of transient flow problems,” , pp. 1–9, 2013.[43] S. T. Zalesak, “Fully multidimensional flux-corrected transport algorithms for fluids,”
Journal of computational physics ,vol. 31, no. 3, pp. 335–362, 1979.[44] M. Rudman, “Volume-tracking methods for interfacial flow calculations,”
International journal for numerical methods influids , vol. 24, no. 7, pp. 671–691, 1997.[45] S. Ii, K. Sugiyama, S. Takeuchi, S. Takagi, Y. Matsumoto, and F. Xiao, “An interface capturing method with a continuousfunction: The thinc method with multi-dimensional reconstruction,”
Journal of Computational Physics , vol. 231, no. 5,pp. 2328–2358, 2012.
46] J. L´opez, J. Hern´andez, P. G´omez, and F. Faura, “A volume of fluid method based on multidimensional advection andspline interface reconstruction,”
Journal of Computational Physics , vol. 195, no. 2, pp. 718–742, 2004.[47] A. Cervone, S. Manservisi, R. Scardovelli, and S. Zaleski, “A geometrical predictor–corrector advection scheme and itsapplication to the volume fraction function,”
Journal of Computational Physics , vol. 228, no. 2, pp. 406–419, 2009.[48] B. Xie and F. Xiao, “Toward efficient and accurate interface capturing on arbitrary hybrid unstructured grids: The thincmethod with quadratic surface representation and gaussian quadrature,”
Journal of Computational Physics , vol. 349, pp.415–440, 2017.[49] R. Scardovelli and S. Zaleski, “Interface reconstruction with least-square fit and split eulerian–lagrangian advection,”
International Journal for Numerical Methods in Fluids , vol. 41, no. 3, pp. 251–274, 2003.[50] W. J. Rider and D. B. Kothe, “Reconstructing volume tracking,”
Journal of computational physics , vol. 141, no. 2, pp.112–152, 1998.[51] M. Sussman and E. G. Puckett, “A coupled level set and volume-of-fluid method for computing 3d and axisymmetricincompressible two-phase flows,”
Journal of computational physics , vol. 162, no. 2, pp. 301–337, 2000.[52] J. Lopez, J. Hernandez, P. Gomez, and F. Faura, “An improved plic-vof method for tracking thin fluid structures inincompressible two-phase flows,”
Journal of Computational Physics , vol. 208, no. 1, pp. 51–74, 2005.[53] R. J. LeVeque, “High-resolution conservative algorithms for advection in incompressible flow,”
SIAM Journal on Numer-ical Analysis , vol. 33, no. 2, pp. 627–665, 1996.[54] L. Jofre, O. Lehmkuhl, J. Castro, and A. Oliva, “A 3-d volume-of-fluid advection method based on cell-vertex velocitiesfor unstructured meshes,”
Computers & Fluids , vol. 94, pp. 14–29, 2014.[55] M. Jemison, E. Loch, M. Sussman, M. Shashkov, M. Arienti, M. Ohta, and Y. Wang, “A coupled level set-moment offluid method for incompressible two-phase flows,”
Journal of Scientific Computing , vol. 54, no. 2-3, pp. 454–491, 2013., vol. 54, no. 2-3, pp. 454–491, 2013.