Latest revision as of 13:39, 26 March 2019

Published in Computers & Fluids, Vol. 141, pp. 90-104, 2016
doi: 10.1016/j.compfluid.2016.04.026

Abstract

In previous works [1,2], the authors have presented a highly efficient extension of the Particle Finite Element Method, called PFEM-2, to solve two-phase flows. The methodology which uses X-IVS [3] to treat convection terms allowing large time-steps was validated for problems where the gravity forces and/or the inertial forces dominate the flow. Although that is the target range of problems to solve with PFEM-2, most of real problems that fall in these categories also includes other flow regimes in certain regions of the domain. Maybe the most common secondary regime is when the surface tension dominates, as an example when drops or bubbles are released from the main flow, and this feature must be taken into account in any complete numerical strategy.

Attending to that, in this work the treatment of the surface tension to PFEM-2 is included. An implicit CSF methodology is employed together with a coupling between the marker function with a Level Set function to obtain a smooth representation of the normal of the interface which allows an accurate curvature calculation. Examples for curvature calculation and isolated bubbles and drops are presented where the accuracy and the computational efficiency are analyzed and contrasted with other numerical methodologies. Finally, a simulation of a jet atomization is analyzed. This case presents the above mentioned features: it is a inertia-dominant flow with a surface tension phenomena on drops and ligaments break up that can not be neglected.

Keywords: PFEM, surface tension, two-phase flows, SCLSVOF

1 Introduction

Solving efficiently multi-phase flows is still an open challenge. Although the dynamics of single phase flows are well understood and can be solved accurately without loss of efficiency, the computational modeling of two or more phases is an underdevelopment field with growing interest. In multi-phase flows the behavior of the fluid at each phase depends on the interface and its shape depends on the flow, then solving this complex coupling is a challenging task.

According to the framework used to derive the formulation, the numerical methods can be split into two main approaches, named Eulerian (fixed framework) and Lagrangian (mobile framework). Former formulations were the first ones to be developed and they provide a natural evolution from single-phase flows since most of Computational Fluid Dynamics (CFD) software are formulated within a fixed framework, while latter formulations offer a more natural choice for simulations in which deformations are not negligible, such as in multi-phase problems.

In the Eulerian strategies, the Finite Element Method (FEM) is a standard tool to solve both structural and fluid problems. In the case of standard FEM, the exact solution cannot be represented in the space spanned by the shape functions, then they will not be able to capture it accurately, but an averaged solution will be obtained. This is particularly important for multi-phase models, since domains composed by different phases usually lead to discontinuities in the properties along the interface, which translates into discontinuities in the unknowns or in their gradients. An alternative to overcome this limitation is to use Enriched Finite Elements [4,5,6,7], which adds degrees of freedom to elements that are cut by the interface in order to capture the part of the solution that escapes from the standard shape function field. Coppola-Owen et al. [8] proposed a simple enrichment functions that is capable of capturing accurately gradient discontinuities (kinks) in the pressure field. Moreover, Ausas et al. [9] proposed a set of three enrichment functions that are able to capture both kinks and jumps in the pressure. Another option in the Eulerian framework is the Finite Volume Method (FVM), which has more followers than the FEM for fluid dynamics. The domain is discretized with cells, and the solution is obtained by calculating fluxes through the faces of each of them. This leads to a formulation that is automatically conservative on the fluxes, unlike FEM.

No matter which Eulerian strategy is used to solve the fluid dynamics, an accurate and efficient simulation of interface evolution is of fundamental importance. For example, in FEM, the use of enriched space is still insufficient to simulate multi-fluids unless it is coupled with a second tool to locate precisely the position of the interface, necessary to build the extra shape functions. It is possible to distinguish two broad classes of computational methods used to describe the evolution of interfaces, namely: interface capturing and interface tracking methods.

Purely Eulerian algorithms, which solve the fluid in a fixed underlying mesh, use capturing methods. In this approach the interface is determined by an implicit function that is advected in the computational domain. Popular methods of this type are the Level Set Method (LSM)[10], which has become widely used when the interface undergoes extreme topological changes, e.g. merging or pinching off; and the Volume of Fluid (VoF) technique[11], which is naturally employed with FVM.

The LSM consists in using a distance function that is convected according to the fluid velocity. This function represents the distance from a point to the interface. By definition, the interface will be located where its value is zero. This level function is variable in the space, but if it has large variations in time, after some time steps it does not represent the distance to the interface anymore, leading to diffusion of the interface and mainly loss of mass. For this reason a reinitialization of the level set must be done to recover a distance function which guarantees that the properties are better conserved. Moreover, an Eulerian advection of the level set function produces large diffusion and requires small time-steps to achieve accurate solutions.

On the other hand, VOF is based on the conservative nature of the FVM, where instead of tracking an interface, it is more natural to save the content of different phases at each cell and define the shape and position from this data later on. The method defines a function that is the fraction occupied by one of the phases in each cell of the domain. Therefore the interface position is not tracked, but the fraction of fluid instead. Once fluids have been convected among cells, the interface position can be reconstructed (accepting some accuracy loss). This exchange between cells, inherited from the conservative nature of the FVM, allows to guarantee mass conservation. This is an important advantage respect to the LSM, in which mass loss is a critical topic which must be addressed and treated. Moreover, the FVM is very robust and is likely to be the most used one in commercial/widespread codes. As an example of application, OpenFOAM®[12] uses this strategy to solve multi-fluid problems.

Formulations clustered in the Lagrangian framework are a more natural choice for simulations where there are large deformations. The original idea, proposed by Monaghan et al. [13] and later works applied to fluid mechanics[14], was a meshless method named Smoothed Particle Hidrodynamics (SPH). Using particles that are advected carrying its own properties over the domain, they are able to almost avoid the numerical diffusion. In the context of incompressible flow, the Lagrangian perspective makes it possible to use a material derivative formulation where the absence of the non-linear convective terms transform the Navier-Stokes system into a transformed linear coupled problem between points and velocities. In the case of multi-phase problems, the calculation of the interface evolution is naturally done using particles [15,16]. However most of Lagrangian formulations have the uncomfortable drawback of requiring a particle position treatment. In the case of meshless methods a constant track of all the moving points must be kept where searching algorithms have to be used to speed up the computational time to calculate the interaction forces. On the other hand, the mesh-based methods must lead with the necessity of constructing or controlling the mesh quality during each time-step of the simulation if the accuracy of the solution has to be maintained. Searching algorithms, evaluation of the mesh distortions or the re-meshing processes are always computationally expensive and it would be interesting to explore the possibility of avoiding those steps.

Alternatives, that combines both Eulerian and Lagrangian tools, have provided to be a good alternative to pure methods. In [17] a pure Eulerian solver for the fluid is used, but Lagrangian marker particles are used to improve the LSM, then the interface tracking. This method proves to be more accurate than the pure Eulerian or pure Lagrangian counterpart in the tracking of the interface. Another option is the named Particle Finite Element Method (PFEM)[18] which consists of using a set of particles that define the nodes of a finite element mesh. Since fluids have no deformation limit, remeshing must be done at each time step. As all Lagrangian methods, the PFEM offers a more natural solution to problems where the particles of the domain can move freely. Unlike LSM, there is no need to recalculate the surface since the location of the interface is obtained trivially; since each particle is associated with a material no extra function is needed. Combining the original idea of Particle in Cell (PIC)[19] where a fixed mesh is used to calculate forces and pressures and moving particles to convect properties, the PFEM method was extended leading to a novel strategy so-called PFEM-2 method[3][20]. This proposal not only tracks material propertiesas density, viscosity, etc., therefore eliminating the need of the non-linear convective term. Also, using an improved explicit integration named X-IVS (eXplicit Integration following the Velocity Streamlines) added to an implicit correction of diffusive terms, there is no limitation in the time step, being the required precision the only bound for the time-step[21]. The enhanced PFEM-2 version to solve multiphase problems, presented in [1] and validated in [2], preserves the large time-step goodnesses of the single-phase strategy, also includes enrichment strategies to capture discontinuities in the pressure gradient, i.e. pressure kinks. However, the range of application of this strategy does not cover an important group of two-phase problems such as those where the surface tension is dominant.

In those problems, a surface tension model must be implemented at the interface being a validated strategy the Continuous Surface Force model (CSF) [22] which is based on an approximation of the interface curvature from the gradient of the marker function. In the case of VoF function, the gradient cannot be calculated accurately since it is a discontinuous step function, and its discrete approximations are known to generate unphysical spurious currents at the interface [23]. Strategies to reduce the spurious currents based on either interface reconstruction or smoothing kernels are available, a literature review can be found in [24], but most of them must be employed only on structured meshes. The coupling achieved by advecting the interface using the conservative VoF function, calculating the interface normal using the smoothed LS function and updating the physical properties from a smoothed Heaviside function is an improved strategy called CLSVOF [25]. In [26], Albadawi et al. present a less expensive option called S-CLSVOF which uses an one-way coupling strategy. This approach was successfully applied on surface tension dominant problems.

The current work proposes a strategy to enlarge the capabilities of PFEM-2 adding the S-CLSVOF method so as to improve the solution of surface tension dominant problems. In Section 1 the governing equations are presented together with the numerical methodology proposed to solve it, doing focus on the interface and surface tension term treatment. Next sections present a set of test properly chosen to show the capabilities of the strategy: starting from pure-convective tests to show the goodness of the Lagrangian framework to transport with neither diffusion nor distortion an arbitrary shape, next surface tension dominant cases where the interface treatment and the method accuracy are quantitatively tested, being PFEM-2 forced to enlarge the time-step where other numerical approaches can not work. Finally, a preliminar simulation of a jet atomization problem is presented. In this type of inertial dominant problem PFEM-2 has demonstrated to be the best numerical option, however a proper treatment of surface tension must be taken into account if an accurate solution of ligaments and droplets formation is searched.

2 Equations and Numerical Formulation

2.1 Governing Equations

The problem to solve is the case of unsteady laminar flow of two immiscible incompressible fluids. Both phases are assumed to be viscous and Newtonian. Isothermal conditions are assumed, and reaction mass transfer and phase transition are not considered. Taking into account all the physical assumptions, both the fluids are then governed by the incompressible Navier-Stokes equation with additional surface tension force along the interface. The governing equations include the continuity, momentum, and interface capturing advection equations which, written in the Eulerian framework, read:

where ${\textstyle \rho }$ is the fluid density, ${\textstyle \mathbf {V} }$ the fluid velocity vector, ${\textstyle \mathbf {g} }$ the gravitational acceleration, and ${\textstyle {\boldsymbol {\sigma }}=-p\mathbf {I} +\mu \left(\nabla \mathbf {V} +\left(\nabla \mathbf {V} \right)^{T}\right)}$ with ${\textstyle p}$ the scalar pressure and ${\textstyle \mu }$ the dynamic viscosity of the fluid. The fluid domain is defined for a single mixture where the function ${\textstyle \lambda }$ is used to distinguish between the two phases. The calculation of the fluid physical properties, the density ${\textstyle \rho }$ and viscosity ${\textstyle \mu }$, vary according to the scalar field ${\textstyle \lambda }$. The surface tension force ${\textstyle \mathbf {F} _{\Gamma }}$ is the concentrated load along the interface, defined as

where ${\textstyle \sigma }$ is the surface tension coefficient between the two phases, ${\textstyle \kappa }$ is the local curvature of the interface, ${\textstyle \delta _{\Gamma }}$ is the Dirac delta function that localizes the surface tension force to point load on the interface and ${\textstyle \mathbf {n} _{\Gamma }}$ as the unit normal to the interface.

Employing the material derivative instead of the temporal derivative plus the convective term, the Equation 1 can be reformulated in its Lagrangian version. In this framework, a kinematic problem has to be solved in order to follow the particle trajectories, which leads to a equation system as:

Due to the multiplication of the velocity and its gradient, the convective term in Equation 1 is non-linear. The presence of this non-linear term demands iterative algorithms to converge to the solution including linearisation techniques to solve the system. On the other hand, the Navier Stokes equations written in a Lagrangian framework lead to a system of linear equations due to the absence of the convective term. Moreover, the resulting system changes the non-symmetric equations in the Eulerian frame into a symmetric and positive definite one.

2.2 Discretization Strategy

The PFEM-2 algorithm consists of solving the equations presented in Section 2.1 using a mixed Eulerian-Lagrangian discretization. This choice leads to a simplification of complex terms, as the convective one, and to a higher accuracy in the results due to lower error in the approximation comparing with its pure Lagrangian or pure Eulerian counterparts [27]. As it was mentioned, since material points move and the configuration changes continuously in time, it is necessary to couple this set of equations with a strategy that solves for the movement of the material points. This is achieved in PFEM-2 by using a set of Lagrangian particles combined with a fixed FEM mesh as is shown in Figure 1. The main advantage of using the Lagrangian particles is that the convection is obtained by simply moving the particles across the space and therefore the system to be solved does not have the convective term. The remaining set of equations will be calculated on the mesh employing a typical fractional step strategy.

Figure 1: Discretization employed in PFEM-2. A cloud of Lagrangian particles are advected over a fixed FEM mesh. Nodal states are projected from neighbor particles states. The neighbor particles are those which are inside the grey region.

It must be remarked that the particles used in the scheme do not represent a fixed amount of mass but rather material points with certain properties and velocity. This allows for different amount of particles to be used depending on the zone to ensure a better accuracy on those areas. Also, it should be noted that in the algorithm presented in this document, the particles are only used to transport the information (solve convective terms of equations). However, in certain cases where the viscosity is low and there is only one fluid phase, it is possible to solve partially the momentum equation in the particles, as explained by [20]. Although this strategy leads to higher accuracy in the cases analyzed in the article, it lacks the generality that is required for the simulation of two-phase problems. Readers interested in a deep explanation of method basis can read [3,20] and the extension to two-phase problems in [1,2] or to fluid-structure interaction in [28].

2.3 Interface Treatment and Surface Tension

In a PFEM-2 simulation, the Lagrangian particles are integrated following a strategy called X-IVAS [3], where the streamlines fixed at time ${\textstyle n}$ are employed to update the particle movements and velocities. In order to extend this approach to track the interface, each particle is initially marked with a sign function ${\textstyle \lambda _{p}}$ depending if it belongs to the first or second phase, and this value is preserved over the particle during the entire simulation guaranteeing boundedness. After X-IVS step, the particle data, i.e. velocity and marker function, must be projected to the nodes to continue with next algorithm steps. Although projection strategies are out of scope of this work (a review and recent improvements of this step can be found in [29]), in order to fix ideas the original algorithm used by PFEM-2 [2] to project a given field ${\textstyle \phi }$ between nodes using subindices ${\textstyle (j)}$ and particles using subindices ${\textstyle (p)}$ is presented in this work, which is of the following form:

where the function ${\textstyle W_{j}}$, associated with the node ${\textstyle j}$, can be either the typical kernel functions used in particle methods such as SPH[13] or the linear shape functions elevated to a power ${\textstyle \alpha }$ (it is ${\textstyle W_{j}(\mathbf {x} )=N_{j}(\mathbf {x} )^{\alpha }}$), while ${\textstyle \mathbf {x} _{p}}$ is the position of the particle ${\textstyle p}$ with state ${\textstyle \phi _{p}}$ and ${\textstyle P}$ is the number of particles in a region around the node ${\textstyle j}$. The region around the node ${\textstyle j}$ can be selected in different ways, being a possibility choosing the zone colored with grey in Figure 1. Mesh nodes thus obtain real values after the projection which can be different to the integer values ${\textstyle \pm {1}}$ that the particles transport. Finally, the interface is defined as the set of points that satisfy the equation ${\textstyle \lambda =0}$.

Once projected over nodes, the function ${\textstyle \lambda }$ has similar properties to a VoF function: the mass is preserved but the discontinuous shape impossibilities an accurate gradient calculation then a poor curvature is estimated which often lead to unphysical flows around the interface when surface tension is included, resulting in unrealistic interface shapes. As it was mentioned, there are several strategies to overcome this limitation, and in this work the approach called S-CLSVOF [26] is selected. An advantage of this approach is that only the ${\textstyle \lambda }$ function is needed to advect (in contrast with CLSVOF), then the initial level set-like function ${\textstyle \phi _{0}}$ is obtained following

The main criterion in choosing this value is to satisfy an initial value of ${\textstyle \phi }$ which is close to the mesh step size. This initial function is a signed function since it has a positive value in the denser fluid and a negative value in the lighter. However, in order to obtain a ${\textstyle |\nabla \phi |=1}$ around the interface, the function is then re-distanced by solving the re-initialization equation:

with the initial condition of ${\textstyle \phi (\mathbf {x} ,0)=\phi _{0}(\mathbf {x} )}$, being ${\textstyle \tau }$ an artificial time discretized with ${\textstyle \Delta \tau =0.1\Delta x}$. Because the re-distancing starts from the initial interface and moving towards both fluids, and we are interested only on the zone around the interface, only few iterations ${\textstyle \phi _{corr}={\frac {\epsilon }{\Delta \tau }}}$ are required, with ${\textstyle \epsilon }$ representing the width around the interface, typically ${\textstyle 1.5\Delta x}$.

After solving 6 the ${\textstyle \phi }$ is now a continuous smooth function around the interface, which helps in determining accurately the interface normal ${\textstyle \mathbf {n} }$ as usual in LSM, it is

Hence, it provides a more precise and smoother interface curvature

Some details are important to remark when Equations 7 and 8 are solved in the FEM framework. An initial strategy is to replace 7 into 8, and to use weighted residuals with the linear piecewise trial functions ${\textstyle \mathbf {N} }$ solving directly for ${\textstyle \kappa }$, it is

However this approach leads to spurious results because ${\textstyle \nabla \phi }$ is discontinuous between the elements. A further option which obtain better results is first obtain a field ${\textstyle \mathbf {\hat {n}} }$ with continuous ${\textstyle \nabla \phi }$ between the elements (linear field in the elements), doing

then, obtain the curvature as usual

with ${\textstyle \mathbf {n} ={\dfrac {\mathbf {\hat {n}} }{|\mathbf {\hat {n}} |}}}$ and ${\textstyle \eta }$ the normal to the boundary ${\textstyle \Gamma }$.

One of the most difficult tasks in front-capturing techniques is to accurately identify the interface to directly impose the term ${\textstyle \mathbf {F} _{\Gamma }}$. This difficulty can be alleviated by interpreting the surface tension as a continuous body force spread across a transition region of thickness avoiding the need of reconstructing the interface explicitly. In this way, the continuum surface force model (CSF) of Brackbill et al. [22] provides an approach to approximate the term of surface tension force ${\textstyle \mathbf {F} _{\Gamma }}$ as a force per unit volume as

where ${\textstyle \delta _{s}}$ is the regularized interface delta function defined as follows

It has been presented [22,30] that if the surface tension term on Equation 12 is discretized explicitly, i.e. the surface tension forces are evaluated on the interface at the previous time step, the stability of the scheme imposes the following restriction on the time step size ${\textstyle \Delta t_{max}}$:

With this restriction the propagation of capillary waves is resolved and their unstable amplification avoided. The Equation 14 can be rather limiting for fine meshes and large surface tension coefficients then is a relevant issue in order to preserve the large time step proposed by PFEM-2. A solution to partially overcome this limitation is to treat the force term 12 implicitly. In this proposal, the surface tension term is included into the implicit calculation of the momentum equation over the mesh using the updated interface position. This helps to extend the time step limitation but is not a fully implicit approach because the interface movement is not coupled with the surface tension imposition. An analysis of the stability of this proposal is presented in Section 4.1.

Finally, the material properties are calculated with the smoothed Heaviside function ${\textstyle H}$

It must be noticed the relevance of the parameter ${\textstyle \epsilon }$ in this strategy. This parameter, which determines the extent of the interface smearing, has been analyzed in other works [26], concluding on the necessity of using the previously mentioned value that preserves a narrow thickness. There are alternative sharp interface methods such as the Ghost Fluid approach [31] which respects jump discontinuities across the interface and avoids an interface thickness. However, in these type of strategies the extension to unstructured meshes is far from straightforward. A FEM framework strategy to treat surface tension without thickness is employing enriched shape functions to treat pressure jumps as proposed by Ausas[9]. In spite of the possibility of capturing jumps, the curvature calculation still being a difficult task. Height functions [32] seems to be the best option, but its formulation for 3d unstructured meshes is still an open challenge. The Laplace-Beltrami formulation [33] appears as an interesting alternative because the curvature does not appear explicitly. However, this strategy also presents drawbacks: it is only accurate with small surface deformations, and requires the computing of the interfacial mesh leading to expensive computations.

2.4 PFEM-2 with surface tension

In order to decouple the unknown fields: velocity and pressure, the projection method known as fractional step [34] is chosen in PFEM-2. This segregated strategy consist on three main steps: velocity predictor, pressure calculation and velocity corrector. The particularity of this predictor step is that the convective term is decoupled from the rest of the momentum equation: the Lagrangian formulation allows to solve it only transporting the particles, which is done employing X-IVS. Then, the particles states are projected to nodal positions, and the remaining terms (including surface tension) are solved implicitly over the mesh, finishing the predictor step with a predictor velocity ${\textstyle {\widehat {\mathbf {V} }}^{n+1}}$ that satisfies the boundary conditions. Pressure calculation ${\textstyle p^{n+1}}$ and velocity correction to obtain ${\textstyle \mathbf {V} ^{n+1}}$ (a divergence free field) are done as usual, but the latter step also includes the particle velocity updating.

It is assumed that all fluid variables are known at time ${\textstyle t^{n}}$ for both the particles and the mesh nodes. Subindexes ${\textstyle ()_{j}}$ y ${\textstyle ()_{p}}$ represent a generic mesh node ${\textstyle j}$ and a generic particle ${\textstyle p}$ respectively (${\textstyle P}$ represents the number of particles). Let ${\textstyle N}$ the finite element linear basis functions. According to this notation, the steps are presented in Algorithm 1, where ${\textstyle \mathbf {x} }$ is a spatial coordinate, ${\textstyle \delta p=p^{n+1}-p^{n}}$, ${\textstyle \delta \mathbf {V} =\mathbf {V} -{\widehat {\widehat {\mathbf {V} }}}}$, ${\textstyle \theta _{p}}$ can be ${\textstyle 0}$ or ${\textstyle 1}$ depending on the pressure restart choice. Also, ${\textstyle \theta _{\mu }}$ can be ${\textstyle 0}$ or ${\textstyle 1}$ depending on the necessity or not of an accurate diffusion calculation when large Fourier numbers are employed, ${\textstyle \mathbf {F} ={\displaystyle \int _{\mathit {\Omega }}N_{j}\left(\rho \mathbf {g} +\mathbf {F} _{\sigma }\right)d{\mathit {\Omega }}}}$, and ${\textstyle \mathbf {M} }$ and ${\textstyle \mathbf {K} }$ are the standard mass and stiffness matrices of any FEM assembling. The computational implementation was done extending the high performing library presented in [21] and each test presented in this work was simulated with that house-made code.

1. Convective Stage:

   ${\displaystyle \left\{\;{\begin{matrix}\mathbf {x} _{p}^{n+1}=\mathbf {x} _{p}^{n}+{\displaystyle \int \limits _{n}^{n+1}}\mathbf {V} ^{n}(\mathbf {x} _{p}^{\tau })d\tau \\{\widehat {\widehat {\mathbf {V} }}}_{p}^{n+1}=\mathbf {V} _{p}^{n}\\\lambda _{p}^{n+1}=\lambda _{p}^{n}\end{matrix}}\;\right.}$

2. Projection Stage:

   ${\displaystyle {\begin{array}{l}\mathbf {M} _{ij}{\widehat {\widehat {\mathbf {V} }}}_{j}^{n+1}=\mathbf {M} _{ip}{\widehat {\widehat {\mathbf {V} }}}_{p}^{n+1}\\\mathbf {M} _{ij}^{L}\lambda _{j}^{n+1}=\mathbf {M} _{ip}\lambda _{p}^{n+1}\end{array}}}$

3. Momentum Stage:

   ${\displaystyle \left(\mathbf {M} ({\dfrac {\rho }{\Delta t}})+\mathbf {K} (\mu )\right){\widehat {\mathbf {V} }}^{n+1}=\mathbf {M} ({\dfrac {\rho }{\Delta t}}){\widehat {\widehat {\mathbf {V} }}}^{n+1}-\theta _{p}\mathbf {G} p^{n}+\mathbf {F} ^{n+1}}$

4. Poisson Stage:

   ${\displaystyle \mathbf {K} \left({\dfrac {\Delta t}{\rho }}+\tau \right)p^{n+1}=\mathbf {B} {\widehat {\mathbf {V} }}^{n+1}+\theta _{p}\left(\mathbf {K} ({\dfrac {\Delta t}{\rho }})p^{n}+\mathbf {B} (\tau ){\boldsymbol {\pi }}^{n}\right)}$

5. Correction Stage:

   ${\displaystyle {\begin{array}{l}\mathbf {M} (\rho )\mathbf {V} ^{n+1}&=&\mathbf {M} (\rho ){\widehat {\mathbf {V} }}^{n+1}-\Delta t\mathbf {G} \left(p^{n+1}-\theta _{p}p^{n}\right)+\theta _{\mu }\mathbf {K} ({\frac {\mu }{\rho }})\left(\mathbf {V} ^{n+1}-{\widehat {\mathbf {V} }}^{n+1}\right)\\\rho _{p}\mathbf {V} _{p}^{n+1}&=&\rho _{p}{\widehat {\widehat {\mathbf {V} }}}_{p}^{n+1}+\sum _{j}\delta \mathbf {V} _{j}^{n+1}N_{j}(\mathbf {x} _{p}^{n+1})\end{array}}}$

3 Interface evolution Tests

This section will deal with an exhaustive validation of the proposed PFEM-2 method to transport arbitrary shapes with neither interface disturbances or mass loss. It is well known that employing the Lagrangian scheme is relatively easy to solve pure-advective problems as presented in this section, but we consider that it is important for the reader to reach a strong conclusion about the goodness of this framework in contrast with the problems observed with the typical Eulerian schemes. The latter are represented by the suite OpenFOAM®which implements a VoF strategy with interface compression. As will be shown, the larger time step is employed, the more relevant are the differences between the frameworks.

3.1 Rigid Body Rotation of Zalesak's Disk

This test consists in the advection of a region composed of a circle with a slot [35]. If the interface tracking is accurate enough, after several revolutions, the shape must remain identical. The computational domain employed is ${\textstyle \Omega \in R^{2}:[0;100]\times [0;100]}$. The advected region is a circle centered at ${\textstyle (50;75)}$ with a radius of ${\textstyle 15}$ and a slot of width ${\textstyle 5}$ and height ${\textstyle 25}$. The velocity field is a rigid body rotation around the center of the domain with a period of 628 time units:

The grid has ${\textstyle 100}$ points in each direction, conforming a cartesian mesh (in the case of PFEM simulations the mesh was split into 20000 triangles). The Courant number used in simulations is aproximately ${\textstyle CFL=4.5}$. Both the initial field and the solution after two revolutions are shown on Figure 2. In the case of PFEM simulation, approximately five particles by element were used. Most relevant OpenFOAM®settings are: SuperBee as the divergence scheme for the linear term in volume fraction advection equation, MULES as the time integration scheme, the number of alphaSubCycles is 20 (to guarantee interface Courant number less than 0.5) and the interface-compression factor cAlpha is set to 1.

Figure 2: Zalesak's disk results after two full revolutions with 100 grid point per direction and ${\displaystyle CFL=4.5}$. The grey region represents the initial condition.

PFEM evolution shows a good agreement with the expected result (shape preservation). Some small errors, which are more evident when the magnitude of velocity is higher, appear due to approximate a curve with a sequence of straight trajectories. Even though in OpenFOAM®simulation the interface-compression method combined with the advection scheme avoids numerical diffusion, they modify the disk shape excessively, finishing in a poor prediction of the final status.

3.2 Single Vortex Case

While Zalesak's disk test is a good indicator of numerical diffusion in an interface-capturing method, it does not test the ability to preserve small scale structures of the fluid flow. A well known test to evaluate the ability of the method to solve structures of different sizes and their evolution is given by the vortex-in-a-box problem introduced by Puckett et al. [36]. The difficulty of this tests is that requires the solution of an interface stretching problem. The computational domain is ${\textstyle \Omega \in R^{2}:[0;1]\times [0;1]}$, where the interest region is a circle centered at ${\textstyle (0.5;0.75)}$ with a radius of ${\textstyle 0.15}$, advected with a velocity field defined by the stream function

being the velocity components

The grid has ${\textstyle 256}$ points in each direction, and the Courant number used in simulations is aproximately ${\textstyle CFL=4.8}$.

The setting employed for each numerical method in this case is almost equal to the previous test, with the only one difference that in OpenFOAM®the interface-compression factor cAlpha is set to 0.25 to give more stability through relaxing in some level the strong sharpness imposition. Using a larger factor, the simulation turns unstable.

Figure 3: Single vortex test using 256 grid points per direction and ${\displaystyle CFL=4.8}$ (${\displaystyle T=8}$). Grey region represents the initial condition.

The results presented in Figure 3 show, for PFEM-2, good agreement with the expected result (shape preservation) after the cycle. Although the first half of the evolution is well captured by OpenFOAM®, the reconstruction of the original shape is not good enough.

3.3 LeVeque Deformation Case

LeVeque [37] proposed a three dimensional incompressible flow field which combines a deformation in the ${\textstyle x-y}$ plane with one in the ${\textstyle x-z}$ plane. This problem can be considered an extension of the previous case, requiring the correct capturing of the stretching phenomenon in three dimensions.

• The computational domain is ${\textstyle \Omega \in R^{2}:[0;1]\times [0;1]\times [0;1]}$.
• The advected region is a circle centered at ${\textstyle (0.35;0.35;0.35)}$ with a radius of ${\textstyle 0.15}$.
• The velocity field is given by

• The grid has ${\textstyle 50}$ points in each direction.
• Courant number used in simulations is aproximately ${\textstyle CFL=4}$.

Since the flow is reversed for ${\textstyle t>T/2}$, after one period the function must return to its original shape. Figure 4 shows that PFEM-2 successfully recovers almost the initial shape, which is a very complicated task for other numerical strategies [38,39].

(a) t = 0	(b) t = T=2	(c) t = T
Figure 4: Snapshots of 3D deformation field test with PFEM-2. There were used 50 points per direction and ${\displaystyle CFL=4}$. Results were smoothed by post-processing purposes.

4 Surface Tension Test cases

In this section several two-phase incompressible tests are presented. The selected cases are focused in problems where the surface tension play an important role. Therefore, the algorithm 2.4 is exhaustive tested in different situations. The preliminary case analyzes the stability of the surface tension modeling measuring the spurious currents and its dissipation level. Being the advection almost negligible, this case allows to show that the Eulerian parts of the algorithm works as other standard codes. Next case consists in a bubble which rises due to buoyancy force under two different regime, one more rigid where the surface tension is stronger and other more inertial where a skirted shape must be found. An analysis of the parasitic currents and mesh convergence is done for the case where the gravity is neglected. Moreover, the same cases simulated with larger time-steps are stable but with more errors when surface tension increases, in contrast to Eulerian algorithms which tends to turn unstable. The second test is a standing wave dominated by capillarity. Although this problem is not the most indicated to be solved by PFEM-2 due the lack of inertial dominance, the method shows good accuracy even using reasonable large time-steps. Finally, a preliminary simulation of a primary atomization of a liquid jet is done. Being an inertia dominated case, large time-steps can be employed but a proper capture of drops and ligaments depends on the local Weber number.

Herein, the efficient distributed-memory implementation presented in [21] and extended to the free-surface treatment (see Algorithm 2.4) is used to simulate each of next cases presented. Also, the numerical parameter ${\textstyle \theta _{p}}$ is set to ${\textstyle 0}$ in every case so as to restart the pressure at each time step to allow larger time-steps and ${\textstyle 3}$ iterations of steps 4 and 5 are done to improve the global first order [2]. On the other hand, the parameter ${\textstyle \theta _{\mu }}$ is set to ${\textstyle 0}$ except when the Fourier number ${\textstyle Fo=\mu \Delta t/\Delta x}$ is greater than ${\textstyle 10}$ where is set to ${\textstyle 1}$. The latter allows to increase the accuracy of the fractional step strategy for highly diffusive problems.

4.1 Stability Analysis

The most critical numerical artifact introduced by the modeling of the surface tension is the generation of spurious currents which appear in the form of vortices around the interface. The employed method, CSF, is not excepted from that drawback. These flows, also named as parasitic currents, are generated solely due to numerical artifacts through the discrete approximation of the interface which acts as a perturbation on the physically smooth interface.

If the surface tension term is discretized explicitly, i.e. the surface tension forces are evaluated at the interface at the previous time step, the stability of the scheme places a stability condition on a time step [22] as

which results in a limiting for fine meshes and large surface tension coefficients. The implicit treatment of surface tension terms is shown to alleviate this restriction [40]. Instead of evaluating the surface tension with the interface at the previous time-step ${\textstyle n}$, the PFEM-2 employs the interface at time ${\textstyle n+1}$. However, the interface movement is not coupled with the surface tension force calculation leading to a sort of semi-implicit scheme. Then, to evaluate the range of stability of this methodology, an analysis similar to that presented by Deshpande et al. [41] is done here.

In order to include the effect of viscosity in the case of low Reynold number, Galusinski and Vigneaux [42] have revisited the time step constrain leading to the following generalized time step criterion

with ${\textstyle \tau _{\mu }}$ and ${\textstyle \tau _{\rho }}$ two independent time-scales depending on the viscosity and density respectively. The constants ${\textstyle C_{1}}$ and ${\textstyle C_{2}}$ are independent of fluid properties and are only solver specific and in this work are determined experimentally from the simulations.

Following [41], the proposed case has a domain of ${\textstyle 1\times 1\mathrm {[mm]} }$ discretized with a uniform grid of ${\textstyle \Delta x=10\mathrm {[\mu m]} }$. Centered is a droplet of radius ${\textstyle R=0.25\mathrm {[mm]} }$. Both density ${\textstyle \rho }$ and viscosity ${\textstyle \mu }$ are the same for the fluid inside or outside the drop and their value depend on the time scale considered. The coefficient of surface tension (${\textstyle \sigma =0.01\mathrm {[N/m]} }$) and simulation time step ( ${\textstyle \Delta t=10^{-4}\mathrm {[s]} }$) are always maintained constant. Gravity is neglected in all simulations. The final simulation time is set to ${\textstyle T_{f}=10\mathrm {[s]} }$.

Figure 5: Stability chart for integration times ${\displaystyle T_{f}=10\mathrm {[s]} }$. Dashed lines represent the boundary between the stable and unstable computations found by Deshpande [41] and by the current work with PFEM-2.

The set of cases simulated covers the values of ${\textstyle \tau _{\rho }/\Delta t}$ and ${\textstyle \tau _{\mu }/\Delta t}$ desired varying ${\textstyle \rho }$ and ${\textstyle \mu }$ properly. The Figure 5 presents the stability charts for the behavior of each test. Three categories of simulations are taken into account:

• Stable: kinetic energy calculated over the entire domain decays and center of mass of the droplet remains fixed.
• Unstable type 1: simulation ends, but the center of mass of the droplet finishes displaced more than the size of one element of the mesh and/or the kinetic energy does not decay.
• Unstable type 2: simulation crashes after completion.

From the results, the constants ${\textstyle C_{1}=0.1}$ and ${\textstyle C_{2}=1}$ can be obtained. It should be noted that this set of simulations has an ideal number of Reynolds of ${\textstyle Re=0}$ which is in the opposite side of a proper application range of Lagrangian strategies. Then, which is actually tested is the performance of the fractional step method employed by PFEM-2 to couple velocity and pressure plus a constant projection/interpolation from/to particles. Comparing with the reference work of Deshpande et al. [41], PFEM-2 is more robust against the density time scale than FVM+VoF, but it is weaker regarding to viscosity time scale. The former is expectable because the density is related to the unsteady and inertia terms. The latter can be understood due to either excessive noise at interface because of the use of particles, which requires more viscosity to dissipate that phenomena than Eulerian strategies and due to failures of the first order fractional step where the iterations do not recompose the solution. Anyway, it has been shown by Deshpande that the generation of spurious currents is only secondary for a moving interface, therefore the time step analysis shown here represents a conservative estimate. It is expected when the convection takes part in the simulation, the advantages of PFEM-2 will be clearer.

4.2 Rising Bubble Case

A widely used surface-tension benchmark is the case of an air bubble rising in a liquid column. Beyond qualitative results, as the bubble shape, Hysing et. al. [43] have presented a set of quantitative results obtained with several CFD multiphase codes solving two cases varying some physical properties. The first one considers a bubble in the ellipsoidal regime which undergoes moderate shape deformation, while in the second one the bubble belongs to the skirted regime and experiences much larger deformation. Both fluids are Newtonian, incompressible and isothermal, with properties listed in Table 1.

In the Fig. 6 the case configuration and the boundary conditions are presented. The initial condition is null velocity with the phase marker imposed as shown. In comparison with other reported works, for example [44], where the initial condition had to be relaxed in order to smooth the interface between the two regions, with the current strategy this pre-processing is not necessary because the initial marker field is imposed over particles then projected to the nodes, obtaining a naturally smoothed field over the mesh which diminishes, but not remove, the typical parasitic current of staircase profiles.

Figure 6: Rising bubble case configuration and boundary conditions.