A coupled volume-of-fluid/level-set method for simulation of two-phase flows in unstructured meshes

A coupled volume-of-fluid/level-set method for simulation of two-phase flows on unstructured meshes N´estor Balc´azara,∗, Oriol Lehmkuhla,b, Llu´ıs Jofrea, Joaquim Rigolaa, Assensi Olivaa,∗ aHeat and Mass Transfer Technological Center (CTTC), Universitat Polit`ecnica de Catalunya - BarcelonaTech (UPC) ETSEIAT, Colom 11, 08222 Terrassa (Barcelona), Spain. E-mail: [email protected] bTermo Fluids, S.L., Avda Jacquard 97 1-E, 08222 Terrassa (Barcelona), Spain. E-mail: [email protected]om Abstract This paper presents a methodology for simulation of two-phase flows with surface tension in the framework of unstructured meshes, which combines volume-of-fluid with level-set methods. While the volume-of-fluid transport relies on a robust and accurate polyhedral library for interface advection, surface tension force is calculated by using a level-set function reconstructed by means of a geometrical procedure. Moreover the solution of the fluid flow equations is performed through the fractional step method, using a finite-volume discretization on a collocated grid arrangement. The numerical method is validated against twoand three-dimensional test cases well established in the literature. Conservation properties of this method are shown to be excellent, while geometrical accuracy remains satisfactory even for the ∗Corresponding author. Fax: +(34) 93 739 89 20 Email addresses: [email protected] (N´estor Balc´azar ), [email protected] (Assensi Oliva ) most complex flows. Keywords: unstructured meshes; level-set method; volume-of-fluid method; finite-volume method; two-phase flow 1. Introduction Numerical simulation of two-phase flows is vital for many engineering and scientific applications, such as combustion, bubbly flow, boiling heat transfer, unit operations in chemical engineering, cooling of nuclear reactors among others. The accurate modeling of interfacial flows is challenging because of the discontinuity in material properties (e.g. density, viscosity), the necessity to account for surface tension force, and due to the fact that the geometry of the interface is not known a priory. For this kind of flows it is critical a precise computation of interfacial quantities such as curvature and normal, which are used to evaluate the surface tension. Errors in the calculated surface tension force will induce non-physical velocities, commonly known as spurious or parasitic currents [32], which can grow with time and so significantly degrade simulation results. Moreover, most of industrial applications are characterized by complex domains, therefore the use of unstructured meshes is advantageous. In order to solve the aforementioned issues many numerical methods have been developed in the past decades. For instance: the front tracking (FT) method [47, 46], level set (LS) methods [30, 42, 29, 1], volume-of-fluid (VOF) methods [18, 48, 24], and hybrid VOF/LS methods [43, 49, 41, 28]. In these 3 methods, two-phase flow i s t reated a s a s ingle fl ow wi th th e de nsity and viscosity varying smoothly across the moving interface which is captured in an Eulerian framework (VOF, LS, CLSVOF, VOSET) or in a Lagrangian framework (FT). Although the idea behind these methods is similar, their numerical implementation may differ greatly. A review of advantages and disadvantages of these techniques in the context of simulation of multiphase flows w ith s harp i nterfaces i s g iven i n [ 48]. I n t he f ront-tracking method [47, 46], a stationary Eulerian grid is used for the fluid flow and the interface is tracked explicitly by a separate Lagrangian grid. This method is extremely accurate but also rather complex to implement due to the fact that dynamic re-meshing of the Lagrangian interface mesh is required [10]. Contrary to LS and VOF method, automatic merging of interfaces does not occur, and difficulties arise when multiple interfaces interact with each other as in coalescence and break-up. In the VOF method [18, 48, 24], the interface is given implicitly by a color function, defined to be the fraction of volume within each cell of one of the fluids. In order to advect the VOF function, the interface needs to be reconstructed using a geometric technique [24]. An advantage of VOF method is the fact that accurate algorithms can be used to advect the interface (e.g. [24]), so that the mass is conserved, while still maintaining a sharp representation of the interfaces [43]. However a disadvantage of the VOF method is the fact that it is difficult to compute accurate curvatures from the volume fraction function used to represent the interface, because it presents a step discontinuity. In level-set (LS) methods [30, 42] the interface 4 is represented by the zero-contour of a signed distance function. The evolution of this function in space and time is governed by an advection equation, combined with a special re-distancing algorithm. One of the advantage of the LS approach is the fact that the interface curvature can be accurately computed, while a disadvantage of this method is that the discrete solution of transport equations leads to numerical error in mass conservation of the fluid-phases. Recently, [29] has introduced a conservative level-set method (CLS) where mass conservation problem is greatly reduced, while an hyperbolic tangent function is employed as the level-set function. Moreover, this approach has been generalized to unstructured meshes by [1, 2, 3], in the framework of finite-volume discretizations. On the basis of the advantages and disadvantages of VOF and LS methods, it can be concluded that they are complementary, so it is an inevitable trend to develop new methods combining VOF and LS approaches. For instance, [43] presented a coupled level-set/volume-of-fluid (CLSVOF) method for computing 3D and axisymmetric incompressible two-phase flows. In the CLSVOF method the curvature was obtained via finite differences of the level set function which in turn is derived from the level set function and volumeof-fluid function. [50] present an adaptive coupled level-set/volume-of-fluid (ACLSVOF) method for interfacial flow simulations on two-dimensional unstructured triangular grids. Another CLSVOF method was implemented by [49] for the numerical simulations of interfacial flows in ship hydrodynamics, where the level set function is re-distanced based on the reconstructed inter5 face with a geometric algorithm, whereas the interface jump conditions were handled by means of a ghost fluid methodology. This method was employed to simulate a gas bubble rising in a viscous liquid and a water drop impact onto a deep water pool. [41] have also presented a coupled volume-of-fluid and level set (VOSET) method, where a distance function is reconstructed from an iterative algorithm and interface is advected by the VOF method. The VOSET method was validated by performing two-dimensional simulations of rising bubbles and dam-break problems. Despite those efforts, to the best of the author’s knowledge most of the aforementioned coupled VOF/LS methods have been designed for regular cartesian meshes, so that their easiness of implementation, capability and accuracy on irregular unstructured meshes is still to be proven. Therefore, the present work is aimed at making progress in the direction of developing an accurate and robust coupled VOF/LS method for simulation of incompressible two-phase flows on twoand three dimensional unstructured meshes, including surface tension effects. Thus, unstructured meshes can be adapted to complex domains, enabling us an efficient mesh distribution in regions where interface resolution has to be maximized, which is in general hard to achieve on structured grids. In the present coupled VOF/LS method, an accurate VOF-PLIC method introduced by [24] is used to advect the interface, while the interface curvature and normals employed to evaluate the surface tension force are computed by using a level-set function. This LS function is reconstructed through a geometrical procedure, based on the computation 6 of the minimum distances between cell centroids and the plane segments provided by the PLIC-VOF method introduced by [24], while the surface tension force is computed in the framework of the continuous surface force model introduced by [5]. Regarding the fluid flow, a classical fractional step method [7] is used to solve the incompressible Navier-Stokes equations, which are coupled with the VOF and LS functions. The Navier-Stokes equations have been discretized by means of the finite-volume method on a collocated unstructured grid arrangement, according to the work introduced in [1]. Numerical results are contrasted against numerical and experimental data from the literature. The outline of this paper is as follows: A summary of the governing equations and numerical methods is given in Section 2. In Section 3 numerical experiments are presented in order to validate the coupled VOF/LS method implemented in this work. These numerical experiments include the simulation of the static droplet test case, twoand three-dimensional buoyant bubbles, co-axial coalescence of two bubbles, and deformation of a drop under shear flow. Finally, the conclusions are presented in Section 4. 2. Governing equations and discretization 2.1. Incompressible two-phase flow The conservation of momentum and mass of two immiscible incompressible and Newtonian fluids is described by the Navier-Stokes equations defined 7 on a spatial domain Ω with boundary ∂Ω: ∂ ∂t(ρkvk) + ∇ · (ρkvkvk) = ∇ · Sk+ρkgin Ωk(1) Sk=−pkI+µk∇vk+ (∇vk)T(2) ∇ · vk= 0 in Ωk(3) Here, Ω = Ω1∪Ω2∪Γ, k={1,2}denote the subdomains associated with the two different fluid phases, Γ = ∂Ω1∩∂Ω2is the fluid interface, ρand µ denote the density and dynamic viscosity of the fluids, vis the velocity field, gis the gravity acceleration, pis the pressure, Sis the stress tensor and Iis the identity tensor. Assuming no mass transfer between the fluids yields a continuous velocity condition at the interface: v1=v2in Γ (4) The jump in normal stresses along the fluid interface is balanced by the surface tension. Neglecting the variations of the surface tension coefficient σgives the following boundary condition for momentum conservation at the 8 interface: (S1−S2)·n=σκnin Γ (5) where nis the unit normal vector outward to ∂Ω1and κis the interface curvature. Eqs. 1-3 and Eqs. 4-5 can be combined into a set of equations for a single fluid in Ω, with a singular source term for the surface tension force at the interface Γ [31, 5, 8]: ∂ ∂t(ρv) + ∇ · (ρvv) = −∇p+∇ · µ∇v+ (∇v)T+ρg+σκnδΓ(6) ∇ · v= 0 (7) where vand pdenote the fluid velocity field and pressure, ρis the fluid density, µis the dynamic viscosity, gis the gravitational acceleration the super-index Trepresents the transpose operator, δΓis a Dirac delta function concentrated at the interface Γ, σis the surface tension coefficient, κis the curvature of the interface and ndenotes the unit normal vector on the interface. Physical properties change discontinuously across the interface: ρ=ρ1H1+ρ2(1 −H1) (8) µ=µ1H1+µ2(1 −H1) 9 with ρ1,ρ2and µ1,µ2the densities and viscosities of the first and second fluids, respectively, whereas H1is the Heaviside step function that is one at fluid 1 and zero elsewhere. In the context of the present VOF/LS method, a volume averaged indicator function will be used in place of H1, as is defined in Eq. 10 and Section 2.7. 2.2. Volume-of-fluid method In the volume-of-fluid method an indicator function fis used to track the interface, f(x, t) =      1 if x∈Ω1 0 if x∈Ω2 (9) with Ω1and Ω2the sub-domains occupied by the fluid 1 and 2 respectively. Discretely, the information effectively stored at the cell ΩPis the volumeaveraged indicator function, namely the volume fraction: fP=RΩPf(x, t)dV RΩPdV (10) where Vis the volume of the cell ΩP. The advection equation for fis given by: ∂f ∂t +v· ∇f= 0 (11) where vis the fluid velocity. 10 2.6. Assessment of the distance function construction The accuracy of the distance function constructed by the algorithm introduced in Section 2.4 is measured. The distance function error E2is defined as E2= 1 ncells ncells X P=1 dP−dexact P2!1/2 (20) where ncells is the number of cells contained in a flagged zone of 3hwidth around the interface, and Pthe index of the current cell ΩP. This test case has been carried out to perform the construction of the signed distance function around a cylindrical bubble of diameter db= 0.5, centered in a unit square domain 2db×2dbdivided by triangular cells, as is illustrated in Fig. 4a. In addition, Fig. 4b shows d(x, t) around three bubbles, with db= 0.4, centered in a unit square domain 2db×2dband triangular mesh with h=db/20. Fig. 5 shows the error E2as function of the grid size h, with p= 1.22 as the convergence order for the signed distance function, around the circle shown in Fig. 4a. Thus, from the present assessment it is demonstrated that the geometrical procedure introduced in Section 2.4 is robust enough to compute d(x, t) with high accuracy, even in presence of multiple interfaces. 2.7. Surface tension force and regularization of physical properties Implementing surface tension in a numerical scheme involves two issues: the curvature κneeds to be accurately calculated, and the resulting pres17 Figure 4: Interface position (black line) and constructed signed distance function (colour lines). (a) Cylindrical bubble of diameter dbon a triangular mesh with h=db/25. (b) Three cylindrical bubbles of diameter dbon a triangular mesh with h=db/20. h E2 0.005 0.025 0.045 10-6 10-5 10-4 10-3 second order first order present (triangular mesh) Figure 5: Error (E2) in distance function computation around a circle of diameter db= 0.5 centered in a unit square domain (triangular mesh). For present results, E2=Chp+ O(hp+1), with p= 1.22 the order of convergence of the algorithm illustrated in Section 2.4. 18 sure jump must be applied appropriately to the fluids. Since the present discretization will be based on the finite-volume integration of the NavierStokes Equations, Eq. (6), the aforementioned problems can be conveniently addressed through the continuous surface force model (CSF) introduced by [5]. Thus, the singular term, σκnδΓ, is converted to a volume force as follows σκnδΓ=σκ(d)∇f(21) where nδΓhas been approximated by the gradient of the volume-averaged indicator function introduced in Eq. 10. From the level-set function, d(x, t), the curvature κis obtained as the divergence of the interface unit normal n: κ=−∇ · n(22) n=∇d ||∇d||(23) To obtain a cell averaged value, the curvature is integrated over each finite volume ΩP: κP=−1 VPZΩP ∇ · ndV (24) 19 Applying the Gauss theorem yields κP=−1 VPZSP n·dA(25) where VPand SPare the volume and surface of ΩPrespectively, while Ais the area vector on SP. Following the work of [25, 1], a wide and symmetric stencil is necessary for accurate evaluation of interface normals in the framework of unstructured meshes, therefore ∇dis computed using the least-squares method with vertex based stencil. The reader is referred to [1] for further details on the application of the least-squares method for gradient evaluation on unstructured meshes. In the context of the present VOF/LS method, physical properties in Eq. 8 are regularized using a volume averaged indicator function, therefore density and viscosity fields are calculated as ρ=ρ1f+ρ2(1 −f) (26) µ=µ1f+µ2(1 −f) with fthe VOF function introduced in Eq. 10. 2.8. Solution procedure for the coupled VOF/LS method Following the work of [1], the Navier-Stokes equations, Eq. 6, have been discretized by means of the finite-volume method on a collocated unstructured grid. A central difference scheme is used to approximate the convective 20 term of momentum equation, Eq. 6, unless otherwise stated; while diffusive terms are centrally differenced. A distance-weighted linear interpolation is used to find the cell face values of physical properties and interface normals, while gradients are computed at cell centroids by using the least-squares method [1]. The velocity-pressure coupling is solved by means of a classical fractional step projection method [7]. The solution procedure used in this work is summarized as follows: 1. Initialize v(xP,0), f(xP,0), d(xP,0), physical properties and interface geometric properties. 2. The time increment ∆t, which is limited by the CFL conditions and the stability condition for the capillary force [5], is calculated by ∆t=C∆tmin h ||v||,ρh2 µ,h ||g||1/2 , h3/2ρ1+ρ2 4πσ 1/2!(27) where C∆t= 0.05 for the current VOF/LS method. 3. The interface is advected by using the PLIC-VOF method introduced in [24] (see Section 2.2). 4. The signed distance function d(x, t) is calculated by the geometric algorithm given in Section 2.4. 5. The curvature is computed by Eq. 22. 6. Physical properties (ρ, µ) are updated by Eq. 26. 21 7. An intermediate velocity v∗is evaluated by ρv∗−ρvn ∆t=−3 2Ah(ρvn)+1 2Ah(ρvn−1)+Dh(vn)+ρg+σκ∇h(φ) (28) where ∇hrepresents the gradient operator, Dh(v) = ∇h·µ∇hv+∇T hv represents the diffusion operator, and Ah(ρv) = ∇h·(ρvv) is the convective operator. For the temporal discretization, explicit AdamsBashforth scheme has been used. 8. The pressure field pis computed by the Poisson equation ∇h·1 ρ∇h(pn+1)=1 ∆t∇h·(v∗) (29) Discretization of Eq. 29 leads to a linear system, which is solved by using a preconditioned conjugate gradient method. 9. The resulting velocity v∗from Eq. (28), does not satisfy the continuity Eq. (7). Therefore it is corrected by vn+1 =v∗−∆t ρ∇h(pn+1) (30) 10. In order to avoid pressure-velocity decoupling when the pressure projection is made on collocated meshes [40, 12], a cell face velocity vf is calculated so that ∇h·v= 0 (see Eq. 7) at each control volume. 22 Namely in discretized form: vf=X q∈{P,F} 1 2vn+1 q+∆t ρ(φn q)(∇hpn+1)q−∆t ρf (∇hpn+1)f(31) where Pand Fare denoting the adjacent cell nodes to the face f. 11. Repeat steps 2-10 until time step required. The numerical algorithms explained in this work have been implemented in the framework of a parallel C++ code called TermoFluids [26]. The reader is referred to [1] for technical details of the spatial and temporal discretizations of the Navier-Stokes equations on collocated unstructured grids. 3. Numerical experiments 3.1. Static drop The first test case is the verification of the stationary Laplace solution for a circular drop with diameter dd. In the absence of viscous, gravitational or external forces, the circular interface with surface tension should remain at rest with the pressure jump at the interface exactly balancing the surface tension force (Laplace’s law): ∆Pexact =σκexact (32) where, the exact curvature is given by κexact = 2/ddfor a circular drop. The correct solution is a zero velocity field and a pressure field that rises from 23 Figure 6: (a) Time evolution of spurious velocities. (b) Pressure profile for different grids, and pressure distribution on Ω for h= 1/100. a constant value of pout =p0outside the drop to a value pin =p0+ 2σ/dd inside the drop. However, at a discretized level, the accurate calculation of the curvature and the balance between the surface tension and pressure jump are not trivial problems, and, as a result spurious currents arise. The computational domain is a square having side lengths of 2ddunits, where dd= 0.5 is the diameter of a bubble positioned at the center of the domain. The coefficient of surface tension, and the viscosity inside and also outside the bubble were all set unity while the densities were given a magnitude of 104. This corresponds to a Laplace number La =ddσρµ−2= 5000, which has been also used by [20]. Present test cases are solved on unstructured meshes of triangular element type, as is summarized in Table 1. For the sake of comparison, the numerical jump in pressure is evaluated 24 Mesh name Number of cells Grid size (h) Cell geometry M11.44 ×1031/25 triangular M25.81 ×1031/50 triangular M32.28 ×1041/100 triangular Table 1: Mesh parameters used in two-dimensional static droplet. σκ(d)∇f σκexact∇f Mesh h L1(v)E(∆p)L1(v)E(∆p) M11/25 1.30 ×10−40.02433 1.12 ×10−40.02508 M21/50 3.19 ×10−50.00651 2.79 ×10−50.00804 M31/100 8.82 ×10−60.00215 7.14 ×10−60.00176 p≈1.94 1.75 1.99 1.92 Table 2: Errors and convergence order (p) for the dimensionless velocity and pressure using the VOF/LS method. Here, {L1(v), E(∆p)}=Chp+O(hp+1) with pthe order of convergence. as follows: E(∆p) = |pin −pout −2σ/dd| 2σ/dd (33) where pin is the pressure inside the drop which corresponds to the maximun pressure on Ω, and pout is the outside pressure which corresponds to the minimum pressure on Ω. Moreover, in order to measure the error in velocity, the following L1error norm is used: L1(v) = 1 Ncells Ncells X k (vk·vk)1/2µ σ(34) which is computed on the whole of the spatial domain Ω. Table 2 shows the errors E(∆p) and L1(v) for different grid sizes (h). As with all Eulerian interface tracking methods there are spurious currents 25 present, however, it can be observed from Fig. 6a that the parasitic currents measured by L1(v) are quite small in magnitude compared with other methods reported in the literature [51, 20], moreover its magnitude tends to a steady state as the time advances. In addition, Fig. 6a and Table 2 show a comparison of L1(v) calculated by the present surface tension model, σκ(d)∇f, against the same model with an exact curvature. It is observed a very slight difference between these results, which confirms the accuracy of the method used to calculate the curvature. Regarding the pressure jump, Fig. 6b and Table 2 illustrate how well the computed pressure fulfilled the Young-Laplace law, Eq. 32, furthermore E(∆p) decreased very rapidly with mesh refinement as is illustrated in Fig. 6b. Finally, data reported in Table 2 have been adjusted by the least-squares method to the function {L1(v), E(∆p)}=Chpto estimate the order of convergence p. The aforementioned results confirm that surface tension model was implemented correctly and it produces accurate results. 3.2. Two-dimensional rising bubble This test case has been solved by [19, 20] in order to determine quantitative reference solutions for the buoyancy-driven motion of a two-dimensional bubble rising in an initially quiescent liquid. The computational setup is illustrated in Fig. 7, where a cylindrical bubble of diameter db= 0.5 is centered in the lower half of a rectangular domain Ω = [0,2db]×[0,4db]. Non-slip boundary condition is applied at the top and bottom boundaries, and free 26 Mesh Grid size (h)Re εRe rsphericity (ζ)εζ r M1d/15 7.04104 0.341% 0.8160 1.18% M2d/20 7.04088 0.339% 0.8125 0.74% M3d/30 7.02213 0.017% 0.8099 0.42% M4d/40 7.01711 −0.8065 − M5d/30 7.02801 −0.8085 − Table 7: Grid convergence and effect of the cell geometry, for Eo = 116, M= 41.1, ηρ= 100 and ηµ= 100. Here εRe r=|Reh=d/40 −Re|/Re,εζ r=|ζh=d/40 −ζ|/ζ. Cell geometry is detailed in Table 6. is defined as ζ=πd2 RΩ||∇f||dV (38) 3.3.1. Grid convergence and cell geometry Fig. 10 and Table 7 show the effect of mesh refinement on the convergence of terminal Reynolds number and sphericity, for the case Eo = 116, M= 41.1, ηρ= 100 and ηµ= 100. There is a very slight difference in results obtained with meshes M3and M4, which correspond to h=d/30 and h=d/40 respectively, therefore mesh M3will be used in discussion of numerical results unless otherwise stated. Moreover, it is observed that estimated errors illustrated in Table 7, are reduced with grid refinement. Regarding the influence of cell geometry on numerical results, Fig. Fig. 10 and Table 7 show a very close agreement in the calculated Reynolds number and sphericity, using the meshes M3and M5, which are formed by triangular prisms and hexahedral volumes respectively. 33 Figure 10: Grid convergence, Eo = 116, M= 41.1, ηµ= 100 and ηρ= 100. Mesh description (M1...M5) in Table 6. (a) Reynolds number. (b) Sphericity (bubble shape). 3.3.2. Effect of convective schemes Numerical simulations have been performed in order to study the influence of the convective scheme used to discretize momentum Eq. 6, on the terminal Reynolds number and bubble sphericity. Following the work of [1], the finite-volume discretization of the convective term of Eq. 6 is based on the use of flux limiters [44], L(θ), defined as L(θ)≡                        1 Central difference limiter (CD), max{0, min{2θ, 1}, min{2, θ}} TVD Superbee limiter, θ+|θ| 1+|θ|TVD Van-Leer limiter, 0 First-order upwind limiter. (39) 34 Figure 11: Effect of the convective scheme used to discretize momentum Eq. 6. Parameters Eo = 116, M= 41.1, ηµ= 100 and ηρ= 100. Mesh M3with h=d/30, details in Table 6. (a) Reynolds number. (b) Sphericity (bubble shape). where θis a monitor variable defined as the upwind ratio of consecutive gradients of the velocity components. The reader is referred to [1] for technical details on the application of flux limiters to discretize the convective term on unstructured grids. Regarding the numerical results, Fig. 11 and Table 8 show that the use of different flux limiters lead to similar results for terminal Reynolds number and sphericity, therefore a CD limiter will be used on the discussion of numerical results unless otherwise stated. However, a TVD-Superbee limiter is also employed in order to avoid numerical instabilities, for instance, in flows that include high Reynolds numbers, ηρ⩾103and ηµ⩾103, or simulations with topology changes. 35 Flux limiter Re sphericity (ζ) CD 7.02213 0.8099 Superbee 7.04084 0.8086 Van Leer 7.03323 0.8089 Upwind 7.01109 0.8139 Table 8: Effect of the convective scheme used to discretize momentum Eq. 6, on the terminal Reynolds number and sphericity, for Eo = 116, M= 41.1, ηρ= 100 and ηµ= 100. Mesh M3with h=d/30. 3.3.3. Bubble shapes and terminal velocity For the sake of comparison, experimental results found in [4, 21, 9] and numerical results reported by [22] are used as reference. Fig. 12a shows the effect of Morton number on the bubble dynamics, given Eo = 116, ηρ= 100, ηµ= 100 and Morton number varying from 5.51 up to 848. It is observed that Re tends to a steady state value for all M, however, as the Morton number increases the overshoot on Re is more pronounced, indicating that the bubble motion has a tendency to reduce its stability. Particularly, the characteristic overshoot of the instantaneous Reynolds number after the bubbles start to ascend is well represented in the cases with M > 5.51. Regarding the mass conservation error of the bubble phase, Fig. 12b proves the satisfaction of this requirement, where a maximum error of O(10−5) is observed. Here, the instantaneous mass is calculated and compared with the initial mass, then mass conservation error is calculated by the expression Mr=|M(t)− M(0)|/M(0) with M(t) = RΩfdV . Fig. 13 illustrates a qualitative comparison of the calculated bubble shapes against experimental images reported by [4], for a wide range of Eo and M, with ηρ= 100 and ηµ= 100. Moreover, Table 9 shows a quanti36 Figure 12: Eo = 116, 1.31 ≤M≤848, ηµ= 100 and ηρ= 100. (a) Time evolution of the Reynolds number. (b) Mass conservation error (Mr=|M(t)−M(0)|/M(0)). Here M(t) = RΩf(x, t)dV . The maximun mass conservation error is O(10−5). tative comparison of terminal Reynolds numbers calculated by the present VOF/LS method against experimental data of [4] and numerical results reported by [22]. In addition Table 10 presents a comparison of VOF/LS results against experimental data taken from the Grace diagram [9], for Eo = 10, 10−3< M < 10, ηρ= 100 and ηµ= 100. Further details of the bubble dynamics, including time evolution of Reynolds number, sphericity, bubble shapes and wake patterns are illustrated in Figs. 14-19. Given the aforementioned results, it is observed that present simulations are in close agreement with previous results from the literature, moreover, numerical stability and accuracy of present VOF/LS method has been proved for a wide range of dimensionless parameters. As further validation for high density and viscosity ratios, a rising bubble 37 Figure 13: Terminal bubbles shapes reported in experiments of [4] (top rows) and computations performed by the present VOF/LS method (bottom rows). Simulations were performed using ηρ= 100 and ηµ= 100. Re Eo M [4] [22] Present Mesh 116 848 2.47 2.317 2.33 M3 116 266 3.57 3.621 3.65 M3 116 41.1 7.16 7.0 7.02 M3 116 5.51 13.3 13.17 13.06 M3 116 1.31 20.4 19.88 19.65 M4 32.2 8.2×10−455.3 52.96 52.84 M3 243 266 7.77 8.397 7.84 M3 339 43.1 18.3 17.91 17.64 M4 Table 9: Present computations for Eo = 116, 1.31 ≤M≤848, ηρ= 100 and ηµ= 100, compared against experimental results from [4] and numerical results from [22]. 38 Re Eo M [9] Present Mesh 10 10−323.6 23.53 M3 10 10−211.7 11.37 M3 10 10−14.9 4.92 M3 10 1 1.7 1.95 M3 10 10 0.6 0.70 M3 Table 10: Present computations for Eo = 10, 10−3≤M≤10, ηρ= 100 and ηµ= 100, compared against experimental results taken from the Grace diagram [9]. Figure 14: Buoyant bubble for Eo = 10, M= 1 ×10−2,ηρ= 100 and ηµ= 100. Mesh M3. (a) Reynolds number and bubble shape evolution. (b) Sphericity and streamlines 39 Figure 15: Buoyant bubble for Eo = 116, M= 266, ηρ= 100 and ηµ= 100. Mesh M3. (a) Reynolds number and bubble shape evolution. (b) Sphericity and streamlines Figure 16: Buoyant bubble for Eo = 116, M= 5.51, ηρ= 100 and ηµ= 100. Mesh M3. (a) Reynolds number and bubble shape evolution. (b) Sphericity and streamlines 40 Figure 17: Buoyant bubble for Eo = 32.2, M= 8.2×10−4,ηρ= 100 and ηµ= 100. Mesh M3. (a) Reynolds number and bubble shape evolution. (b) Sphericity and streamlines Figure 18: Buoyant bubble for Eo = 243, M= 266, ηρ= 100 and ηµ= 100. Mesh M3. (a) Reynolds number and bubble shape evolution. (b) Sphericity and streamlines 41 Figure 19: Buoyant bubble for Eo = 339, M= 43.1, ηρ= 100 and ηµ= 100. Mesh M4. (a) Reynolds number and bubble shape evolution. (b) Sphericity and streamlines Figure 20: Buoyant bubble for Eo = 39.4, M= 0.065, ηρ= 714 and ηµ= 6670. Mesh M3. (a) Dimensionless velocity, vb·ey(gd)−1/2. (b) Sphericity and streamlines 42 Figure 25: (a) Computational setup and initial condition. (b) Droplet deformation parameter, D. (c) Mesh configuration, triangular prism cells. coalescence process are illustrated in Fig. 24c. As the bubbles rise, a liquid jet is formed behind the leading bubble, which induces a severe deformation in vertical direction of the following bubble. Then, once the two bubbles are approaching, the trailing bubble accelerates because the suction by the top bubble. As time progresses, the two bubbles start to touch, leaving a mushroom-like structure. Finally, the thin liquid film between bubbles is squeezed out and ruptured, completing the coalescence process. Numerical predictions match fairly well in terms of bubble shapes with experimental results reported by [6]. Moreover, present results are consistent with VOF simulation from [48] and level-set simulations reported by [1]. 3.5. Deformation of a droplet in a shear flow A spherical drop of diameter dis located at the center of a computational domain x∈[0,8], y∈[0,4], z=∈[0,8], without effect of gravity force, 49 Figure 26: (a) Vorticity contours (ex· ∇ × v) on the droplet surface and on the plane z−yat x= 0. (b) Time evolution of droplet shapes on the plane z−yat x= 0. Here ηρ= 1, ηµ= 1, Re = 0.1, 0.05 ⩽Ca ⩽0.3. Red line used for steady state shape. Figure 27: (a) Droplet sphericity versus time. (b) Capillary number Ca versus Taylor deformation parameter D. Present VOF/LS method (red symbol); BVOF [27]; −theory [45]; 4Lattice-Boltzmann method [54]; 5boundary integral method [52]; 50 as shown in Fig. 25a. The opposite velocities ±Uare imposed on the top and bottom walls, periodic boundary condition is applied in xdirection and Neumann boundary condition in ydirection. The initial condition at time t= 0 is a drop with spherical form and the initial velocity field is linear inside the computational domain. Computations have been performed using an unstructured mesh formed by 1.94 ×106triangular prism cells with grid size h=d/30, as is illustrated in Fig. 25c. This mesh was generated by a constant step extrusion of a two-dimensional triangular grid along the y−axis, being hthe step size. The deformation behaviour of the droplet is determined by the Reynolds number (Re) and the capillary number (Ca). The viscosity is given by the Reynolds number Re ≡ρc˙γd2 4µc (48) where the shear rate is defined by ˙γ= 2U/Lz. The capillary number is given by Ca ≡˙γµdd 2σ(49) This dimensionless parameter is a measure of the relative effect of the shear stress versus the surface tension across the fluid-fluid interface. The Reynolds number used in these test cases is Re = 0.1, while the capillary number is in the range 0.05 ⩽Ca ⩽0.3. The same viscosity and density are specified for 51 both drop fluid and continuous fluid, thus ηρ=ρc/ρd= 1 and ηµ=µc/µd= 1, where the sub-index cis used for the continuous fluid and dthe droplet fluid. The droplet shapes at steady state and some vorticity distributions are illustrated in Fig. 26. The interface becomes ellipsoidal and its deformation and rotation are larger as the capillary number increases. A theoretical solution was derived by [45] to predict small distortions of the droplets from the spherical form at slow speeds, on the hypothesis of Stokes flow. This result show that the droplet is distorted into an ellipse where the deformation parameter given by D= (L−B)(L+B)−1(see Fig. 25b) is linearly changed with the capillary number (Eq. 49). Here Land Bdenote the semi-major and semi-minor axes of the ellipse, as is illustrated in Fig. 25b. Fig. 27a shows the time evolution of droplet sphericity for 0.05 ≤Ca ≤0.3, whereas in Fig. 27b are plotted the Taylor deformation parameter (D) versus capillary number (Ca). Here, it is observed a close agreement between present computations using VOF/LS method, against previous results from the literature [45, 27, 54, 52]. Moreover, for small capillary numbers the shape in steady state is close to the theoretical predictions of [45], while the theory underestimates the droplet deformation parameter for large capillary numbers, as shown Fig. 27b. 52 4. Conclusions A coupled VOF/LS method, which combines the advantages and overcomes the disadvantages of both techniques, has been proposed for computing incompressible two-phase flows on unstructured meshes. From the comparison of the present numerical simulations against experiments and numerical data from the literature, it is possible to conclude that this method is enough robust to perform high accurate computations of interfacial flows with surface tension. Moreover, an error less than O(10−5) in the mass conservation property of the fluid phases is achieved because a VOF-PLIC method is used for interface advection, while an accurate computation of interface curvature and surface tension is performed by means of a level-set function reconstructed from a geometrical algorithm. In addition, numerical stability of the present unstructured VOF/LS solver has been proved for a wide range of dimensionless parameters, including simulations with high density and high viscosity ratios, and interfacial flow with topological changes. Altogether, these validations demonstrate that the present VOF/LS approach and the developed code for simulating two-phase flows on collocated unstructured meshes can be used for practical applications. 5. Acknowledgments This work has been financially supported by the Ministerio de Econom´ıa y Competitividad, Secretar´ıa de Estado de Investigaci´on, Desarrollo e Innovaci´on, Spain (ENE2011-28699), and by Termo Fluids S.L. N´estor Balc´azar 53 acknowledges financial support in form of a doctoral scholarship of the Agencia Espa˜nola de Cooperaci´on Internacional para el Desarrollo (AECID), Spain. Three-dimensional simulations were carried out using computer time provided by PRACE (project 2014112666) through the MareNostrum III supercomputer based in Barcelona, Spain. References [1] Balc´azar, N., Jofre, L., Lehmkhul, O., Castro, J., Rigola, J., 2014. A finite-volume/level-set method for simulating two-phase flows on unstructured grids. International Journal of Multiphase Flow 64, 55-72 [2] Balc´azar, N., Lemhkuhl, O., Rigola, J., Oliva, A., 2015. A multiple marker level-set method for simulation of deformable fluid particles. International Journal of Multiphase Flow 74, 125-142 [3] Balc´azar, N., Lemhkuhl, O., Jofre, L., Oliva, A., 2015. Level-set simulations of buoyancy-driven motion of single and multiple bubbles. International Journal of Heat and Fluid Flow 56, 91-107 [4] Bhaga, D., Weber, M.E., 1981. Bubbles in viscous liquids: shapes, wakes and velocities, J Fluid Mech 105, 61-85 [5] Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A Continuum Method for Modeling Surface Tension, J. Comput. Phys. 100, 335-354. 54 [6] Brereton, G., Korotney, D., Coaxial and Oblique Coalescence of Two Rising Bubbles, in Dynamics of bubbles and Vortices Near a Free Surface, vol. 119, ASME, New York, 1991 [7] Chorin, A.J., Numerical solution of the Navier-Stokes equations. 1968. Math. Comput. 22, 745-762. [8] Chang, Y.C., Hou, T.Y., Merriman, B., Osher, S., A level-set formulation of Eulerian interface capturing methods for incompressible twophase flows. 1996. Journal of Computational Physics 124, 462-488. [9] Clift, R., Grace, J.R., Weber, M.E., Bubbles, Drops and Particles. Academin Press, New York, 1978. [10] Deen, N.G., Van Sint Annaland, M., Kuipers, J.A.M. , 2009. Direct numerical simulation of complex multi-fluid flows using a combined front tracking and immersed boundary method, Chemical Engineering Science 64, 2186-2201. [11] Engquist, B., Tornberg, A.K., Tsai, R., 2005. Discretization of Dirac delta functions in level set methods, Journal of Computational Physics 207, 28-51. [12] Felten, F.N., Lund, T.S., 2006. Kinetic energy conservation issues associated with the collocated mesh scheme for incompressible flow, J. Comput. Phys. 215, 465-484. 55 [13] Gottlieb, S., Shu, C., 1998. Total Variation Dimishing Runge-Kutta Schemes, Mathematics of Computations 67, 73-85. [14] Gueyffier, D., Li, J., Nadim, A., Scardovelli, R., Zaleski, S., 1999. Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows, J. Comput. Phys. 152, 423-456. [15] Haselbacher, A., Vasilyev, O., 2003. Commutative discrete filtering on unstructured grids based on least-squares techniques, J. Comput. Phys. 187, 197-211. [16] Hadamard, J.S., 1911. Mouvement permanent lent dune sphere liquide et visqueuse dans un liquide visqueux, C. R. Acad. Sci. 152, 1735. [17] Harmathy, T.Z., 1960. Velocity of large drops and bubbles in media of infinite or restricted extend, AIChE J. 6, 281-288. [18] Hirt, C., Nichols, B., 1981. Volume of fluid (VOF) method for the dynamics of free boundary, J. Comput. Phys. 39, 201-225 [19] Hysing, S., Turet, S., Kuzmin, D., Parolini, N., Burman, E., Ganesan, S., Tobiska, L., 2009. Quantitative benchmark computations of two-dimensional bubble dynamics, International Journal for Numerical Methods in Fluids 60, 1259-1288. [20] Hysing, S., 2012. Mixed element FEM level set method for numerical simulation of immiscible fluids, J. Comput. Phys. 231, 2449-2465. 56 [21] Hnat, J.G., Buckmaster, J.D., 1976. Spherical cap bubbles and skirt formation, Phys. Fluids 19, 182-194 [22] Hua, J., Stene, J., Lin, P., 2008. Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method, J. Comput. Phys. 227, 3358-3382 [23] Joseph, D., 2003. Rise velocity of a spherical cap bubble, J. Fluid Mech. 488, 213-223 [24] Jofre, L., Lehmkuhl, O., Castro, J., Oliva, A., 2014. A 3-D Volume-ofFluid advection method based on cell-vertex velocities, Computers & Fluids 94, 14-29 [25] Kothe, D.B., Rider, W.J., Mosso, S.J., Brock, J.S., 1996. Volume Tracking of Interfaces Having Surface Tension in Two and Three Dimensions, AIAA 96-0859. 25. [26] Lehmkuhl, O., Perez-Segarra, C.D., Soria, M., Oliva, A., 2007, A new Parallel unstructured CFD code for the simulation of turbulent industrial problems on low cost PC cluster, Proceedings of the Parallel CFD 2007 Conference, pp.1-8. [27] Li, J., Renardy, Y.Y., Renardy, M., 2000. Numerical simulation of breakup of a viscous drop in simple shear flow through a volume-offluid method, Physics of fluids 12, 269-282. 57 [28] L´opez, J., G´omez, P., Hern´andez, J., Faura, F., 2013. A two-grid adaptive volume of fluid approach for dendritic solidification, Computers & Fluids. 86, 326-342. [29] Olsson, E., Kreiss, G., 2005. A conservative level set method for two phase flow, J. Comput. Phys. 210, 225-246. [30] Osher, S., Sethian, J.A., 1988. Fronts propagating with curvaturedependent speed: Algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79, 175-210. [31] Peskin, C.S., 1977. Numerical analysis of blood flow in the heart, J. Comput. Phys. 25, 220-252. [32] Raesi, M., Mostaghimi, J., Bussmann, M., 2010. A volume-of-fluid interfacial flow solver with advected normals, Computers & Fluids 39, 1401-1410. [33] Renardy, Y., Renardy, M., 2002. PROST: A Parabolic Reconstruction of Surface Tension for the Volume-of-Fluid Method, Journal of Computational Physics 183, 400421. [34] Rodrigue, D., 2001. Generalized correlation for bubble motion, AIChE Journal. 47, 39-44. [35] Ryskin, G., Leal, L.G., 1984. Numerical solution of free-boundary problems in fluid mechanics. Part 2. Buoyancy-driven motion of a gas bubble through a quiescent liquid, J. Fluids Mech. 148, 19-35. 58