Hierarchical x-fem applied to n-phase flow

7th Workshop on Numerical Methods in Applied Science and Engineering (NMASE 08) Vall de N´uria, 9 a 11 de enero de 2008 c LaC`aN, www.lacan-upc.es HIERARCHICAL X–FEM APPLIED TO N–PHASE FLOW Sergio Zlotnik1∗and Pedro D´ıez2 1: Group of Dynamics of the Lithosphere (GDL) Institute of Earth Sciences “Jaume Almera”, CSIC Llu´ıs Sol´e i Sabar´ıs s/n, 08028 Barcelona, Spain 2: Laboratori de C`alcul Num`eric, Departament de Matem`atica Aplicada III Universitat Polit`ecnica de Catalunya Campus Nord UPC, 08034 Barcelona, Spain Palabras clave: multiphase flow; level set methods; enrichment; eXtended Finite Element Method (X–FEM) Resumen. In this work we proposed an extencion of the level set technique to track any number of free surfaces. This extension is based in a hierarchical ordering of several level set functions. To complete the X–FEM approach, the enrichment via partition of the unity method is also extended. The ridge function, base of the enriched interpolation, is restated to include several level sets and the hierarchy between them. 1 INTRODUCTION Despite the term “multi flow” is widely used in the level set community, most works using level sets for tracking free surfaces limit the number of phases to two. In these works, the sign of a level set describes the phases location. There are some exceptions of “multi–phase” or “n–phase” models which can handle n>2. For example, the work of Tan and Zabaras [7] combines level sets with some features of front tracking methods to model the microstructure evolution in the solidification of multi–component alloys. In this work each component is defined by a level set function: the sign limits the solid–liquid interface. Two algorithms to simulate triple junctions where the motion of the interfaces depends on surface tension and bulk energies were proposed by Zhao et. al. [9] and Ruuth [5]. Both use several level sets to track interfaces. They use as many level sets as materials; to prevents overlapping or vacuum, some artificial constraints are added to the model. 2 Problem statement We will focus on unsteady incompressible n–phase (n>2) viscous flows of immiscible fluids which can be described by the Stokes equations in its quasi–static version (inertia 1 Sergio Zlotnik and Pedro D´ıez term neglected). This is a common approach in geophysical modeling, where creeping (very slow) flow arises. The governing equations can be written as ∇·(η∇su)+∇p=ρg,(1a) ∇·u=0 (1b) where uis the velocity, ηthe viscosity, pthe pressure, ρthe density, and gthe gravitational acceleration vector. The symetric gradient operator ∇sis defined as 1/2(∇+∇). Density and viscosity fields are constant on each phase, leading to discontinuities across all interfaces. As the equation (1) is quasi–static, it does not contain any explicit time dependence and the transient character of the solution is due to the motion of the phases. The location of the different phases is described by a collection of level set functions. The level sets represent material properties and they are consequently transported by the motion of the fluid. Thus the evolution of each one of the level sets, describing phase locations, is determined by pure advection equation ˙ φ(i)+u·∇φ(i)=0 (2) where uis the velocity field, solution of the Stoke’s problem (1), and φ(i)is the level set number i. 3 Describing a n–phase fluid with n−1level sets The level set technique is widely used in two and three dimensions to track the interface between materials in two–phase flow problems, see for example [8, 1, 10]. 3.1 Two phases with a single level set The location of the interface between two materials can be described using a level set function φ(1).Thesuperscript (1) denotes the number of level set and it will be useful when 3 or more phases were described. Despite it is not necessary for the discussion in this Section, it is included here to use the same notation as in the following Sections. The sign of the level set φ(1) describes a partition of the simulation domain Ω in two subdomains Ω1and Ω2using the following sign convention φ(1)(x,t)=⎧ ⎨ ⎩ x∈Ω1if φ(1)(x,t)>0 xis on the interface if φ(1)(x,t)=0 x∈Ω2if φ(1)(x,t)<0 (3) where xstandsforapointinΩandtis the time. The interface location is the set of points where the level set field vanishes. An example of partition is shown in Figure 1. Initially the level set φ(1) is defined as a signed distance to the interface. Far enough from the interface, φ(1) is truncated by maximum and minimum cutoff values. The resulting level 2 Sergio Zlotnik and Pedro D´ıez set function describes the position of the interface independently of the computational mesh, thus the same mesh can be used trough the entire simulation avoiding remeshing procedures. Ω Ω1 Ω2 Figure 1: One level set function splits the domain in two subdomains corresponding to the different phases. 3.2 Tracking more than two phases: hierarchy of level sets One level set allows for describing only two phases (two subdomains). To include a third subdomain Ω3a second level set function φ(2) is needed. We propose to assign a hierarchy to the level set functions: the subdomain Ω1is determined by the first level set φ(1) as φ(1)(x,t)=x∈Ω1if φ(1)(x,t)>0 x/∈Ω1if φ(1)(x,t)<0(4) The curve where the level set φ(1)(x,t) equals zero is the interface between the first phase and the rest of the domain. That is, either the second or the third phase. The remaining part in the simulation domain (x∈Ω\Ω1) is split by the second level set φ(2) as for x/∈Ω1,φ (2)(x,t)=x∈Ω2if φ(2)(x,t)>0 x∈Ω3if φ(2)(x,t)<0(5) determining the location of the second and third sub domains. Note that the second level set does not have any influence where the first level set is positive. The first level set is “prior to” —or has upper hierarchy than— the second level set. Figure 2 shows the partition of the domain by two hierarchical level sets into three subdomains. This hierarchy phase description can be extended to the general case of nphases being tracked by n−1 level sets. The level set number i,φ(i), defines the location of the phase ias follows for x/∈ i−1  j=1 Ωj,φ (i)(x,t)=x∈Ωiif φ(i)(x,t)>0 x/∈Ωiif φ(i)(x,t)<0(6) 3 Sergio Zlotnik and Pedro D´ıez Ω Ω1 Ω2Ω3 φ(1) φ(2) Figure 2: Two hierarchical level sets describe three material sub domains. The second level set φ(2) acts only where the first level set φ(1) is negative. Dotted line represents the level set with lower hierarchy eclipsed by the first level set. for all i=1...n−2. The less hierarchical level set φ(n−1) determines the location of the last two phases in the remaining space as for x/∈ n−2  j=1 Ωj,φ (n−1)(x,t)=x∈Ωn−1if φ(n−1)(x,t)>0 x∈Ωnif φ(n−1)(x,t)<0(7) In this approach the positive part of the i–th level set defines the material subdomain Ωiand the negative region have to be partitioned by the level sets with less hierarchy. Figure 3 illustrates a partition into four subdomains by three level sets. Ω Ω1 Ω2 Ω3 Ω4 φ(1) φ(2) φ(3) Figure 3: Three hierarchical level sets allows for describing four material phases. The last level set φ(3) acts only where the first two level set are negative. A shorter description of the domains defined by each hierarchical level set can be done by means of the McCauley brackets φ=1/2(φ+|φ|). 4 Sergio Zlotnik and Pedro D´ıez The domain Ωidescribed by the level set φiis Ωi= suppφi i−1  j=1 −φj. 4 Enrichment with X–FEM The interface described by a level set does not need to conform with mesh edges, that leads to elements with different material properties in its interior. In X–FEM approach, the interpolation in those elements is enriched allowing the expected gradient discontinuity across the interface. The interpolation of velocity uin enriched elements is composed by the standard finite element part, plus an enriched part. The second involves additional degrees of freedom ajand its associated interpolation functions Mj uh(x,t)= j∈N uj(t)Nj(x)+  j∈Nenr aj(t)Mj(x)(8) where Nis the set of standard finite element velocity degrees of freedom and Nenr is the set of enriched degrees of freedom. The Nenr set evolves trough time and needs to be recomputed at each time step after level set movement. The pressure field pis enriched in a similar way. The interpolation function Mjis constructed as the product of standard nodal shape functions and a ridge function R Mj(x)=Nj(x)R(x).(9) The Rfunction is based on the level set and has a “crest” just over the interface between materials. Several different ridge functions have been proposed in the literature (see for example [1, 2]). In the next sections a ridge function based on several hierarchical level sets is proposed. It is based on the ridge for two–phases used by Mo¨es in [2]. 4.1 Two phases, one Ridge In two–phase simulation the only interface is described by the only level set. The enriched elements are those which are crossed by φ(1), and the ridge function can be constructed as R(x)=  j∈Nenr |φ(1) j|Nj(x)− j∈Nenr φ(1) jNj(x).(10) Note that this ridge function vanishes in the element edges not crossed by the level set. Thus, the solution between enriched and non–enriched elements conforms naturally. The plot of such a ridge is shown in Figure 4. 5 Sergio Zlotnik and Pedro D´ıez node 1 Ridge Level set node 2 node 3 element interface φ(1) 2 φ(1) 3 φ(1) 1 Ω1Ω2 Figure 4: Ridge function Rbasedononelevelsetφ(1). 4.2 Ridge function for two level sets When two hierarchical level set are used two things have to be redefined. Firstly, the detection of enriched element has to include the level set hierarchy. Secondly, the triple junction case, where two level sets intersects, has to be take into account. The enriched element detection can be expressed as follows: find the elements crossed by φ(1) and the elements crossed by φ(2) with φ(1) negative. This statement can be easily encoded, in the code repository we provide a highly vectorized MATLAB function named crossedByLevelSet which accept any type element in any number of dimensions and returns if one element has to be enriched. The ridge function Rin enriched elements cross by only the k–th level set is defined as in the previous case. The function r(k)is the ridge associated with level set k r(k)(x)=  j∈Nenr |φ(k) j|Nj(x)− j∈Nenr φ(k) jNj(x)(11) where kis one or two for elements crossed by φ(1) or φ(2), respectively. The ridge function is R=r(k). In the triple junction case, where two interfaces cross simultaneously one element, the ridge function has to take into account both level sets and the hierarchy between them. In this case Ris defined as R(x)=r(1)(x)+r(2)(x)C(1)(x) (12) where the cutoff C(1) function introduces the level set hierarchy C(1)(x)=⎧ ⎨ ⎩ 1ifφ(1)(x,t)≤0 0ifφ(1)(x,t)≥m(1) f(1)(x)otherwise (13) Here m(1) is the minimum positive nodal value of the level set φ(1) in element, and f(1)(x) is defined as f(1)(x)=1−φ(1) m(1) 2.(14) 6 Sergio Zlotnik and Pedro D´ıez The cutoff C(1) function restricts the second ridge r(2) to the region where the first level set φ(1) is negative, and smoothly reduces the value of R(2) where φ(1) is positive. Composite ridge R= r(1) + r(2)C(1) (d) Cutoff function C(1) (c) Ridge r(2) level set 2 φ(2) 3 φ(2) 2 φ(2) 1 (b) node 3 node 2 node 1 Ridge r(1) level set 1 φ(1) 3 φ(1) 2 φ(1) 1 (a) Ω1 Ω2 Ω3 node 3 node 2 node 1 (e) Figure 5: Building of a composite ridge based in two hierarchical level sets. Plots (a) and (b) show the level sets and the simple ridge based on it. The cutoff function C(1) based on the first level set is shown in plot (c). The construction of a ridge function inside an element crossed by two level sets is shown in Figure 5. Note that obtained ridge shown in panel (d) conforms with its three neighbor elements: edge 2–3 is shared with a non enriched element and the ridge is zero. Edge 1–3 is shared with an element crossed only by the first level set and the ridge on this edge takes the same values as r(1). Finally, the ridge in edge 1–2 depends on values of both level sets, in this case the neighbor element is crossed by the same level sets too. 4.2.1 General case The ridge of the last example can be extended to the general case where nlevel sets simultaneously crossing one element. In that case, the ridge function is extended in the following form R(x)=r(1)(x)+ n−1  i=2 r(i)(x)C(i−1)(x) (15) with the C(i)cutoff function C(i)(x)=⎧ ⎨ ⎩ 1ifφ(i)(x,t)≤0 0ifφ(i)(x,t)≥m(i) f(i)(x)otherwise (16) 7 Sergio Zlotnik and Pedro D´ıez the f(i)function f(i)(x)=1−φ(i) m(i)2(17) and the value m(i)=min(|φ(i)|).(18) Note that when the ridge function is constructed with an level set which not crosses the element, is becomes zero. Thus in Equation (15) all ridges can be added together. 5 Numerical examples The proposed n–phase approach is used to simulate some gravitational Raleigh–Taylor instabilities. The model is composed by three fluid immiscible materials. The mechanic of the problem is governed by Equation (1). The driving force in all the presented models is the gravity; the lower buoyant layers are less dense than the overlying layers and a diapir develops. The initial configuration of the following models is composed by three materials located as shown in figure 6. The upper layer is ten times denser than the two lower materials. As the lower materials have different viscosity the formed diapir looses its vertical axis of symmetry. 01 0 0.5 0.5 1 1.5 g Figure 6: Initial configuration composed of three materials. The upper material is denser than the other two. The viscosity of the materials is indicated on figure 7. The evolution of four models with different viscosity contrast between the two lower layers is shown in Figure 7. The (a) row corresponds to a model where all materials have thesameviscosityη= 1. In this conditions the two buoyant materials behave as a unique fluid and a standard symmetric diapir develops. The second row in Figure 7 shows the evolution of the diapir when the viscosity of the right lower material is five times the viscosity of the left material. In this case the symmetry is lost. The evolution of the left half of the model is similar to the (a) row while the right half of the model is controlled by the viscosity contrast between the right material and the overburden layer. The models of 8 Sergio Zlotnik and Pedro D´ıez the third and fourth rows have a viscosity contrast between the two lower layers of 10 and 100, respectively. The very viscous right material of the last model is almost stopped, while the left material develops the diapir alone. The generated flow change its main pattern during evolution. In the early stages (1st and 2nd snapshots) the high viscosity of the right material inhibits the movement in the right half and the flow is concentrated in the left part of the domain. This flux inclines the diapir to the left. Once the material gains enough height to loose the influence of the viscous layer (last two snapshots), the main flow moves to the right half of the model because there are more space facilitating the return flow. This right flux inclines the diapir to the right. 6CONCLUSIONS We have proposed and tested a methodology to extend the level set technique to track any number of free surfaces. The extension is based in a hierarchical ordering of several level set functions. To complete the X–FEM approach, the enrichment via partition of the unity method is also extended. The ridge function, base of the enriched interpolation, is restated to include several level sets and the hierarchy between them. REFERENCES [1] J. Chessa and T. Belytschko. An extended finite element method for two–phase fluids. Transactions of the ASME, pages 10–17, 2003. [2] N. Mo¨es, M. Cloirec, P. Cartaud, and J. F. Remacle. A computational approach to handle complex microstructure geometries. Computer Methods in Applied Mechanics and Engineering, 192:3163–3177, 2003. [3] S. Osher and R. Fedkiw. Level set methods: an overview and some recent results. Journal of Computational Physics, 169:463–502, 2001. [4] S. Osher and J.A. Sethian. Front propagating with curvature dependent speed: algorithms based on hamiltonjacobi formulations. Journal of Computational Physics, 79:12–49, 1988. [5] S. J. Ruuth. A diffusion–generated approach to multiphase motion. Journal of Computational Physics, 145:166–192, 1998. [6] J.A. Sethian and P. Smereka. Level set methods for fluid interfaces. Annual Review of Fluid Mechanics, 35:341–372, 2003. [7] L. Tan and N. Zabaras. A level set simulation of dendritic solidification of multicomponent alloys. Journal of Computational Physics, 221:9–40, 2007. [8] G. J. Wagner, N. Mo¨es, W. K. Liu, and T. Belytschko. The extended finite element method for rigid particles in Stokes flow. International Journal for Numerical Methods in Engineering, 51(3):293–313, 2001. 9