A locally extended finite element method for the simulation of multi-fluid flows using the Particle Level Set method

Accepted Manuscript A locally extended finite element method for the simulation of multi-fluid flows using the Particle Level Set method K. Kamran, R. Rossi, E. O˜ nate PII: S0045-7825(15)00190-5 DOI: http://dx.doi.org/10.1016/j.cma.2015.05.017 Reference: CMA 10630 To appear in: Comput. Methods Appl. Mech. Engrg. Received date: 30 July 2014 Revised date: 5 May 2015 Accepted date: 6 May 2015 Please cite this article as: K. Kamran, R. Rossi, E. O˜ nate, A locally extended finite element method for the simulation of multi-fluid flows using the Particle Level Set method, Comput. Methods Appl. Mech. Engrg. (2015), http://dx.doi.org/10.1016/j.cma.2015.05.017 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 A locally extended finite element method for the simulation of multi-fluid flows using the Particle Level Set method K. Kamrana,∗, R. Rossia,b, E. O˜ natea,b aCentre Internacional de M´etodes Num´erics en Enginyeria (CIMNE), Gran Capit´an s/n, 08034 Barcelona, Spain bUniversitat Polit`ecnica de Catalunya, Barcelona, Spain Abstract The simulation of immiscible two-phase flows on Eulerian meshes requires the use of special techniques to guarantee a sharp definition of the evolving fluid interface. This work describes the combination of two distinct technologies with the goal of improving the accuracy of the target simulations. First of all, a spatial enrichment is employed to improve the approximation properties of the Eulerian mesh. This is done by injecting into the solution space new features to make it able to correctly resolve the solution in the vicinity of the moving interface. Then, the Lagrangian Particle Level Set (PLS) method is employed to keep trace of the evolving solution and to improve the mass conservation properties of the resulting method. While the local enrichment can be understood in the general context of the XFEM, we employ an element-local variant, which allows preserving the matrix graph, and hence highly improving the computational efficiency. Keywords: Two-phase flows, Enrichment, Discontinuous pressure, XFEM, Particle Level Set. 1. Introduction The simulation of multi-fluid problems faces two main challenges. The first one is related to the difficulty in approximating kinks and jumps in the simulated field, and the second one is connected to the difficulty in accurately capturing or tracking the interface between different fluids [1]. If the interface-tracking approach is used, the computational domain adapts itself to the shape and position of the interface. The mixed Lagrangian-Eulerian approach in [2], the space-time approach in [3] or the pure Lagrangian particle finite element (PFEM) approach in [4] are among this family. The way the interface is tracked is on itself a parameter of paramount importance in the treatment of the jump/kink in unknown fields. The main drawback of interface-tracking methods is the need for regenerating the computational mesh at the distorted zones. Recent developments in the PFEM seem to circumvent this problem by convecting particles and then projecting and solving on the fixed mesh [5]. ∗Corresponding author. Tel.: +34 93 401 7399; Fax: +34 93 401 6517 Email addresses: [email protected] (K. Kamran), [email protected] (R. Rossi), [email protected] (E. O˜ nate) The alternative approach to the treatment of the interface, called interface-capturing method, is based on the idea of capturing the position of the interface by keeping trace of some marker field. Typical examples of this approach are the Volume of Fluid method (VOF) and the Level Set method. Interface-capturing methods have gained more popularity as they easily deal with merging and breaking-up of the interface. In the Level Set method, the interface is represented implicitly as the zero-level of a smooth function and typically the signed distance function is the candidate. Since the values of the level-set function are needed exclusively in the vicinity of the interface, it is often observed that this equation is only solved there [6, 7]. Despite many efforts in using high order spatial and temporal schemes to evolve the distance function, it needs to be reinitialized frequently as it soon ceases to be a distance function. Unfortunately, it appears to be difficult to conserve the interface during the reinitialization process which often leads to the breakdown of the mass conservation. Various solutions have been proposed to overcome this problem. In [8, 9] the volume-of-fluid (VOF) method is coupled with the Level Set method to obtain a second-order technique which is generally superior to either method alone. The Particle Level Set (PLS) method introduced in [10] uses Lagrangian marker par- Preprint submitted to Comput. Methods Appl. Mech. Engrg May 4, 2015 *Manuscript Click here to download Manuscript: XFEM.pdf Click here to view linked References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 ticles to rebuild the level set in regions which are underresolved. This is often the case for flows undergoing stretching and tearing. Particles are seeded near the interface and contain a local measure of the interface position that is used to correct the level set function error due to convection or reinitialization. Geometric mass-preserving redistance method proposed in [11, 12] and developed for unstructured meshes, avoid a uniform mass-conserving correction. This method computes a node-wise correction on a narrow band close to the interface considering that mass loss/gain concentrates on the interface zones with higher curvature. The Discontinuous Galerkin (DG) method has also been applied to improve mass conservation of the level set method. The quadrature-free DG method developed in [13] for the conservative form of the level set equation remains stable even if the level set diverges from a signed distance function. Unfortunately, as we anticipated before, even when the position of the interface in space is known, a second challenge arises: accurately approximating the fields of interest within the cut elements close to the interface. Sharp definition of the interface results in cut elements. If jumps in material properties are large or interface forces are involved, the conventional FEM can not capture the possible jump/kink in the pressure/velocity fields. One solution is to assign a numerical thickness to the interface that provides a smooth transition in the material properties at both sides of the interface [14]. This type of methods requires constant interface thickness during time. However, if a sharp interface is desired, cut elements need to be enriched. In [15] pressure is locally enriched at the interface to capture a kink in the pressure field. In [16] two enrichment functions, one for each side of the cut, are introduced at the elements crossed by the interface. These functions are local to each element, linear on each side of the interface, discontinuous along the interface and zero at the element nodes. Similar to [15], these enrichments are condensed before assembly. Beside local enrichment, global one, and in particular the XFEM have been also widely used to model multi-fluid flow [17–20]. Except for the intrinsic XFEM [21], all versions of the XFEM add enrichments to the global system and therefore the graph of the system needs to be updated as the interface moves. In this work we present a local enrichment only for the pressure field and couple it with the PLS method. The basis for the enrichment are taken from the shifted XFEM [22], or its equivalent [23], and presented for triangular and tetrahedral elements. Basically inside each element we add one enrichment for each node and then condense them all at the element level. This enrichment can be seen as an extension to the ones proposed in [15] and [16] that locally add one or two DOFs, respectively. We then couple this enrichment with the PLS method to capture the interface. The overall system is presented in the residual-based variational multiscale stabilized form similar to the one proposed in [15, 20]. The content of this paper is presented as follows. In section 2 the incompressible two-fluid problem is presented. Our proposal for the local enrichment is introduced in Section 3, and Section 4 is devoted to the PLS. Comparison between the XFEM and our proposal and also the effect of the PLS in the results are presented through different numerical examples in Section 5. 2. Problem statement The incompressible Navier-Stokes equations governing the two-phase motion in a domain Ωare written as: ρ(∂tu+u·∇u)−∇·(2µ∇su)+∇p=ρb ∇·u=0.(1) This set of equations is completed by the appropriate Dirichlet and Neumann boundary conditions. The domain Ωis split into two parts, denoted by +and −, with the following material properties: (ρ(x), µ(x)) =((ρ+, µ+) if x∈Ω+ (ρ−, µ−) if x∈Ω−. Note that a jump in density at the interface Γproduces a discontinuity in the pressure gradient, and a jump in viscosity causes a discontinuity in pressure. Surface tension can also produce a jump in the pressure field at the interface. In this work we do not consider the effect of surface tension and therefore we have the balance of the internal forces at the interface as: (σ−−σ+)·n=0. The internal stress σat each domain has the form σ= −pI +2µ∇su. The variational equivalent of (1) is to find (u,p)∈ 2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 V×Qsuch that: ZΩ ρ(∂tu+u·∇u)·vdΩ + ZΩ 2µ∇su:∇svdΩ −ZΩ p∇·vdΩ = ZΩ ρb·vdΩ ZΩ q∇·udΩ = 0, ∀(v,q)∈V×Q. A stabilized variational form can be obtained using the Algebraic Sub-grid Scale (ASGS) method [24] or the Finite Calculus [25] approach. Nonlinearities are best handled in an implicit fashion using a residual-based Newton-Raphson strategy. Temporal discretization is performed by a second order Bossak scheme [27] (or see page 543 in [26]). The semidiscrete stabilized variational form of this problem is: Ru=ZΩGu·vhdΩ + ZΩ 2µ∇sun+1 h:∇vhdΩ −ZΩ pn+1 h∇·vhdΩ +X e τ1ZΩe un h·∇vh·(Gu+∇pn+1 h)dΩe +X e τ2ZΩe∇·un+1 h∇·vhdΩe=0 Rp=ZΩ qh∇·un+1 hdΩ +X e τ1ZΩe (Gu+∇pn+1 h)·∇qhdΩe=0.(2) In (2) the term Guis given by: Gu=ρ(∂tun+1 h+un+1 h·∇un+1 h−bn+1), and the stabilization parameters are, τ1=h2 e 4µ+2ρhe|ue|and τ2=µ+0.5ρhe|ue|, where heand ueare the characteristic elemental length and the velocity vector, respectively. Accurate solution of (2) requires a modification in the finite element space, to capture the possible kink or jump expected at the solution. 3. Enrichment When the interface cut elements on a fixed mesh, the jump in density causes a kink in the pressure field and hence it induces a jump in its gradient at the interface Figure 1: Two-fluid hydrostatic flow with a jump in density. Linear elements can not represent a kink in the pressure field inside an element. A) The interface passes through the elements, b) the jump in density in elements cut by the interface, and c) the exact solution (dotted line) and the FE solution (solid line). position (Figure 1). It is clear that for simple elements, triangles in 2D and tetrahedra in 3D, the linear approximation of pressure can not represent the expected kink in pressure. The pressure field does not belong anymore to the standard one, and therefore one way to capture it is to add enrichments to the standard field in these elements. Remark.The method described in this work is mainly applied to physical problems dominated by gravitational forces, i.e. small errors in the pressure can lead to large errors in the velocity. Therefore, an improvement of the pressure approximation is more beneficial than modifying the velocity approximation. Numerical studies performed in [19] suggest that it is not advisable to enrich the velocity approximation space as it does not improve the results significantly. Severe convergence problems and an increase in the required number of iterations have been also reported as the consequences of the velocity enrichment. We define the enriched pressure field at each cut element as: pe h= Nnode X i Ne ipe i+ Nnode X i Mi(pe enr)i.(3) Here (pe enr)iare the elemental pressure enrichments associated to the corresponding enrichment basis Mi. These basis are defined at each node of the cut element as: Mi(x,t)=Ni(x,t)·[ψ(x,t)−ψ(xi,t)],(4) with ψ(x,t) being the global enrichment function and Ni(x,t) is the standard FE shape function for node i. This definition for Mi(x,t) has the property that Mi(x,t) is zero in the mesh nodes and takes its maximum at the interface. Zero support at the nodes facilitates interpreting the new added DOF as local to the cut element rather than global. This definition for Mi(x,t) ensures 3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Kronecker-δproperty of the overall approximation. It means that at each given node i, pressure is only interpolated with its own nodal value because any other node jhas zero support at node i. The enrichments do not affect explicitly the nodal approximation of pressure as they are zero at the nodes of the mesh. These basis are known as the XFEM shifted basis [19] in the literature. Two definitions are proposed for ψin the literature. In case of strong discontinuities (jump in pressure) sign enrichment is defined as: ψsign(x,t)=sign(φ(x,t)) =         −1 if φ(x,t)<0 0 if φ(x,t)=0 1 if φ(x,t)>0. (5) where φ(x,t) is the level set function that is defined as the signed distance function whose zero matches the interface position. Its exact definition is presented in Section 4. Similarly for the weak discontinuities (kink in pressure fields) Mo¨ es et al. [28] and Hansbo & Hansbo [23] separately proposed an abs-enrichment shape function that is defined as: ψabs(x,t)= Nnode X i|φi|Ni(x,t)−| Nnode X i φiNi(x,t)|. Figure 2 shows the enrichment shape functions, Mi, for weak and strong discontinuities in a triangular element. Remark.Numerical results in [19] suggest that even when no surface tension is considered, sign-enrichment provides better results than abs-enrichment although the latter seems sufficient. Using the sign-enrichment, a larger approximation space for the pressure is provided which also allows capturing strong discontinuities in pressure. On the other hand, we know that a discontinuity in pressure not only appears due to the surface tension effect, but also due to viscosity jumps [29]. Therefore, In this work we only consider the signenrichments: ψ(x,t)=ψsign(x,t). As already mentioned, these enrichments have zero support at the nodes and therefore can be interpreted as local to each cut element. The matrix form of the semidiscrete stabilized form (2), using equations (5), (4) and (3), at each cut element is written as:        A B BTD       U Penr != F 0!.(6) Uis the vector of nodal velocities and pressures of the (a) (b) 0 0.5 (c) ψabs 1 -1 (d) ψsign -1 0 (e) Msign 1 -1 0 (f) Msign 2 0 1 (g) Msign 3 0 1 (h) Mabs 1 0 1 (i) Mabs 2 0 1 (j) Mabs 3 Figure 2: A) triangular element cut by the interface. B) sub-elements used for the numerical integration. ×shows the integration points. C) ψabs, d) ψsign, e-g) sign-enrichment basis to capture a strong discontinuity and h-j) abs-enrichment basis to capture a weak discontinuity. element, U={u,p}, and Penr contains the elemental enrichments for pressure. Following this arrangement, matrix Acontains all terms related to the original DOFs, U, matrix Bcontains terms related to the original DOFs and enrichments, and matrix Dcontains terms related only to the enrichments, Penr. Note that, to be consistent, the enrichment is also applied to the test function of pressure. Vector Frepresents the body force term. Exact computation of the integrals in (2) requires a modification in the quadrature rule. To this end, each triangular element (tetrahedral in 3D) is split into sub-elements and for each sub-element the same integration rule as for the non-cut elements is used (see Figure 2b). Now we proceed to the condensation by computing first the enriched pressures Penr from the second equa- 4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 tion in (6) as: Penr =−D−1BTU, and then plug it into the first equation in (6) to obtain, (A−BD−1BT)U=F, at each cut element. The condensation process is valid as long as the matrix Dis invertible. Recalling equation (2), Dis the Laplacian of the enriched pressures computed as: D=τ1ZΩe∇penr ·∇qenr dΩe. From the definition of Min (4), it is clear that, no matter where the interface is, the ∇Miis always bounded and actually constant in each side of the interface. As the result, Dis always invertible unless if for some Mithe support tends to zero. Remark.The enrichment shape functions with small support typically occur when the interface is close to a node or edge and lead to an ill-conditioned system matrix. Two approaches to treat this problem are given in [30–32]. In this work we use the blocking criteria proposed in [30] in which for each cut element enrichment is considered if the ratio between the minimum and maximum sub-volumes is bounded by a userdefined constant C  1, i.e. min(Ve 1,Ve 2) max(Ve 1,Ve 2)<C.(7) Note that doing condensation, the cut element ends up having the same DOFs, U, as the rest of the elements and therefore the graph of the matrix remains unchanged as the interface evolves. This is not the case for the conventional XFEM where the enrichments are global. Figure 3 compares our proposed local enrichment and the XFEM. In local case, Figure 3a, only the cut elements are engaged and the additional DOFs are represented with black circles inside each element. On the other hand in the XFEM enrichment scheme, Figure 3b, three types of subdomains are usually distinguished. The first one, Ωenr, is composed of elements cut by the interface which have all nodes enriched. The second subdomain, Ωpenr, has all elements that have at least one node enriched, and the last subdomain is the collection of elements with standard DOFs and without enrichment. In XFEM, as the interface moves, new elements are cut and so new global DOFs need to be considered. This means that the graph of the matrix assembly need to be recomputed, which is quite costly for large problems. Remark.Note that the XFEM provides a continuous pressure field at the nodes and at the inter-element boundaries, while the local enrichments are only continuous at the nodes (they are zero at the nodes) and not at the inter-element boundaries. This irregularity in pressure is still acceptable because pressure, as appears in (2), belongs to the L2(Ω) space and does not require control over the derivatives. (a) Elemental enrichment (b) XFEM Figure 3: Schematic view of the condensed elemental enrichment and the XFEM. A) local enrichment. Solid points represent the condensed pressure DOFs at each cut element. B) subdomains containing enriched, partially enriched and standard (non-enriched) elements as well as the enriched nodes (solid nodes) in the XFEM terminology. Remark.Similar local enrichments that add one enrichment for the cut element or two, one at each side of the interface, already exist in the literature. One of such modifications [15] suggested to interpolate the pressure as: pe h= Nnode X i Ne ipe i+Ne enr pe enr, where Nnode is the number of nodes per element and Nenr is the new enrichment function added just for the 5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 elements cut by the interface. This function has the constant gradient at each side of the interface and has zero support at the nodes. Nenr is defined by the level set values φat the element nodes as: Nenr =0.5(−| Nnode X i Ne iφi|+ Nnode X i Ne i|φi|). This definition for Nenr is quite similar to the definition for ψabs shown in Figure 2c. Note that here this function is directly used as the enrichment shape function while in XFEM it is used to define Mi. The enrichment proposed in [15] suffers from instabilities. The idea proposed here of capturing the kink in the pressure field with just one added DOF perfectly works for 1D problems. However, for 2D and 3D situations one enrichment is not sufficient to capture the correct pressure field and actually one enrichment per cut edge seems more appropriate. Intuitively, as long as the edges are cut at the same height, one enrichment is sufficient but when deviation occurs and therefore each edge is cut at a different height, more than one enrichment is needed. As for 3D problems the number of cut edges can reach up to four in each tetrahedra, a generic enrichment that covers all cases would be an enrichment at each node. Our proposed local XFEM satisfies this requirement as all nodes are enriched, and at the same time takes advantage of the condensation to avoid change in the graph of the matrix as is common in conventional XFEM. In the following section we review an interface capturing technique based on the Level Set method. No iteration is considered between the multi-fluid solver and the interface capturing one. At each time step, first the interface is moved and then one step of the predictor multi-corrector fluid solver is performed. This is in contrast to some more sophisticated coupling strategies such as those presented in [33] and [19]. 4. Particle Level Set (PLS) method The underlying idea behind the Level Set method is to represent an interface as the zero-level set of a higher dimensional function φ(x,t). This function is scalar and substantially reduces the complexity of describing the interface, especially when undergoing topological changes such as pinching and merging. The level set function φ(x,t) is defined to be a smooth function that is positive in one region and negative in the other. The motion of the interface is determined by a velocity field, u, which can be a function of position, time and geometry of the interface, or be given externally, for instance as the solution of the Navier-Stokes equations (1). The advection equation for the interface evolution is: φt+u·∇φ=0.(8) This level set equation is a first order hyperbolic PDE and only needs to be solved near the interface. The most common choice for the level set function φis the signed distance to the interface so that |∇φ|= 1. This ensures that the level set is a smoothly varying function well suited for high order accurate numerical methods. There are several techniques to solve the level set equation (8) in space and time [6, 7]. Despite the high order temporal and spatial approximations of the level set equation, instabilities may appear when the level set cease to be a signed distance function. This situation occurs at the presence of large topological changes at the interface vicinity, which are quite common in practical problems. One solution is to reshape the level set function to a distance function. This method, called reinitialization, has been shown to stabilize the numerical instabilities and is performed frequently during the evolution of the interface. Reinitialization is performed at each time step by solving to steady state (as fictitious time τ→ ∞) the equation: φτ+sgn(φ0)(k ∇φk −1) =0.(9) where sgn(φ0) is a one-dimensional smeared out signum function [14]. Unfortunately, one of the major drawbacks of the reinitialization process is the difficulty in preserving the original location of the interface, often leading to breakdown in the conservation of mass. To overcome this problem of mass loss with the level set method (Figure 4a), various solutions have been proposed [10, 11, 13, 34]. The PLS [10] method uses Lagrangian marker particles to rebuild the level set in regions which are underresolved. This is often the case for flows undergoing stretching and tearing. Two sets of massless marker particles are placed near the interface with one set, the positive particles, in the φ > 0 region and the other set, the negative particles in the φ < 0 region (Figure 4b). It is unnecessary to place particles far from the interface and this greatly reduces the number of particles needed in a simulation. The region near the interface could be considered as the region covered by all elements that have at least one corner with the distance inferior to three times the element size. In this work 10 particles are seeded 6 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 in each tetrahedral element, which is different from the 4d(d is the dimension) particles proposed in the original work for the cubic elements [10]. Considering that each cubic element could be divided into 6 tetrahedra the choice of 10 particles is reasonable. The placement of the massless particles around the zero-level set for this test can be seen in Figure 4b. In Figure 4c the PLS solution after one revolution is shown. (a) (b) (c) Figure 4: Particle Level Set (PLS) method [10]. A) the mass gain/loss of the standard level set solution after one revolution. B) placement of massless positive (blue) and negative (red) particles around the interface for the Zalesak test. C) PLS solution after one revolution. For the purpose of interface reconstruction a radius rpis assigned to each particle as the function of its distance to the interface. This radius is bounded by minimum and maximum values based upon the mesh size near the interface (0.1d<rp<0.5d,dis the mean mesh size). Note that this radius provides a local measure to the interface position. Once particles are seeded and their radius adjusted, the level set equation (8) is integrated forward in time. Then, particles are advected with the evolution equation dxp dt =u(xp),(10) where xpis the position of the particle and u(xp) is its velocity. The particle velocities are interpolated from the velocities on the underlying grid. Remark.In this work we use the second order Crank- Nicolson method for the time stepping of the level set equation (8). The particle evolution equation (10) is integrated in time using sub-stepping technique in which 10 sub-steps are considered. In smooth well resolved regions of the flow where the level set method is highly accurate, the particles do not drift a non-trivial distance across the interface allowing us to maintain a high order accurate level set solution. Instead, in under-resolved regions we find particles that are on the wrong side of the interface by more than their radius. We define a particle as escaped particle only if it crosses the interface by more than its radius. Escaped particles are used to reconstruct the level set function in under-resolved regions. To this end, a local level set function is defined for each escaped particle by means of the radius associated to the particle as: φp(x)=sp(rp− k x−xpk),(11) where spis the sign of the particle, i.e. ±1. These level sets are only defined locally on the corners of the element containing the escaped particle and can be seen as the particle predictions of the values of the level set function on the corners of the element. Any variation of φfrom φpindicates possible errors in the level set solution. The escaped positive particles are used to rebuild the φ > 0 region and the escaped negative particles to rebuild the φ≤0 region. For example, take the φ > 0 region and an escaped positive particle. Using equation (11), the φpvalues of the grid points on the element containing the particle are calculated. Each φpis compared to the local values of φand the maximum of these two values is taken as φ+. This is done for all escaped positive particles creating a reduced error representation of the φ > 0 region. That is, given a level set φand a set of escaped positive particles E+, we initialize φ+with φ on all grid points and then calculate φ+=max ∀p∈E+(φp, φ+). Note that this operation is done only for the escaped positive particles and φp’s are only calculated on the elements containing the escaped particles. Similarly for the negative region, φ≤0, we initialize φ−with φand then calculate φ−=min ∀p∈E−(φp, φ−). φ+and φ−will not agree due to the errors in both the particle and the level set solutions as well as interpolation errors, etc. We merge φ+and φ−back into a single level set value by setting φto the value of φ+or φ−which is least in magnitude at each grid point, φ=(φ+if |φ+|⩽|φ−| φ−if |φ+|>|φ−|. The minimum magnitude is used to reconstruct the interface (instead of, for example, taking an average), since it gives priority to values that are closer to the interface. The PLS method is also capable of correcting the errors due to the reinitialization. During the reinitializa- 7 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 tion step, particles are not moved to conserve the zerolevel set and then any error in the reinitialization scheme is corrected by them. After the reinitialization and correction steps, the radii of the particles are adjusted according to the current value of φ(xp). In summary the order of operation in one time step of the PLS method is as follows: •Solve the level set equation (8) and move the particles. •Correct errors in the level set function using particles. •Reinitialize the distance function and again correct errors using particles. •Adjust the particle radii. •Reseed particles, if necessary, to have a uniform distribution. Our proposed local XFEM is loosely coupled with the PLS method to perform numerical examples as shown in the following section. At each time step first the interface is moved, and then one step of the predictor multi-corrector fluid solver is performed. 5. Numerical examples All the following examples are solved on tetrahedral meshes and no surface tension is considered. Regarding the PLS method, the maximum seeding distance is limited to 3h, where his the maximum edge size. The constant used in the enrichment criteria of Eq. (7) is taken as C=10−4. 5.1. Rayleigh-Taylor instability A single mode Rayleigh-Taylor instability is considered to study the interface capturing technique coupled with the enrichment method. This instability occurs when a heavy fluid is accelerated into a light one due to gravity forces. Various authors have used XFEM or other stabilized methods coupled with the Level Set method and hexahedral elements for the simulation of this phenomenon [20, 35]. A rectangular domain Ω = [L×H] with L=0.5mand H=4mis extruded in the third direction by one element thickness h=0.015625 m. The domain is uniformly discretized by 49152 tetrahedra (32 ×256 ×1 cubes, each of which divided into 6 tetrahedra). The instability is triggered by a cos form deviation of the amplitude 0.05 at the horizontal interface between two fluids, z=0.05cos(2πx). The top heavy fluid has a density ρ+=3kg/m3and a dynamic viscosity µ+= 0.0135 kg/(ms) and the bottom light fluid has ρ−= 1kg/m3and µ−=0.0045 kg/(ms). The gravitational acceleration of magnitude g=−10 m/s2in the zdirection is considered. Similar dimensionless Reynold number, Re =√gHLρ+/µ+=√gHLρ−/µ−=500 and Atwood number, A=ρ+−ρ− ρ++ρ−=0.5, as in [20, 35] are obtained. Slip boundary conditions are considered on the side walls, no-slip condition on the bottom wall and zero pressure is prescribed at the upper wall. No surface tension is considered and a fixed time step of ∆t=0.01s is used. Figure 5 shows the evolution of the interface at different instances along the formation of the instability. Our results are compared well with the results obtained in [20, 35] with hexahedral meshes of the same size. The maximum mass fluctuation is approximately 0.12%. In [20] a mass fluctuation of 0.10% and in [35] a mass fluctuation of 0.07% was observed for a secondorder level-set scheme. (a) 0.5 s (b) 0.8 s (c) 1.1 s (d) 1.25 s (e) 1.5 s Figure 5: Rayleigh-Taylor instability. A-c) shows the evolution of the interface until it breaks. D-e) shows the particles involved in the interface reconstruction concentrated near the interface. Note that the mesh has low resolution to capture the underlying structure. 5.2. Sloshing tank The ability of the proposed enriched method to model the free surface behavior is studied in this example. 8