A new incompressibility discretization for a hybrid particle MAC grid representation with surface tension
Introduction
Incompressible free surface flows are of great interest due to their wide scope of applications ranging from small scale droplet dynamics to large scale dam breaking. Simulation of incompressible free surface flows is a challenging problem because the location of the free surface and hence the shape of the fluid domain is time dependent. Surface tension forces exacerbate the situation by imposing a strict time step restriction. A number of approaches have been proposed to alleviate this restriction. [24] derived an implicit formulation of surface tension by looking at the changes of curvature over time. This was simplified in [26] via the use of the Laplace–Beltrami operator. This simplified approach was also implemented for finite volume models in [41]. The approaches of [49], [61], [64] split the Laplace–Beltrami operator into a standard Laplacian term and an extra term allowing for a semi-implicit approach treating the standard Laplacian implicitly. There are also approaches that improve stability of VOF-based surface tension schemes [54], [35], [51]. A main difficulty faced by these approaches is the lack of an explicit representation of the interface when deriving an implicit model for surface tension. [43] discretized the surface tension force on a Lagrangian mesh attempting to avoid this difficulty. They coupled the Lagrangian surface tension force with the pressure on the Eulerian grid using an interpolation matrix. This significantly improved the stability, but due to the indirection of the interpolation matrix the surface tension force does not achieve an exact balance with the incompressible pressure and thus still suffers from instability with larger time steps.
Motivated by [43], a Lagrangian description of the interface is desirable when deriving an implicit model for surface tension because interface positions, curvature and force derivatives are readily accessible. This motivates the use of a front tracking method [59], [56]. But front tracking alone does not completely alleviate the issues because one still has to balance surface tension with internal pressure forces, and this coupling is still difficult to address when the pressure forces are described on an Eulerian grid. An alternative solution is to take a fully Lagrangian approach to the problem and solve the incompressibility condition directly on particles. These particle-based methods have gained much attention for incompressible flow problems, see e.g. [14], [13], [25], [34], [62]. However, [55] pointed out that these solvers are expensive due to the inefficiency of solving the Poisson equation on an irregular particle mesh. They also pointed out that mapping particles to a background grid and solving for incompressibility on that grid is more efficient as in [23], [10], [36]. However, the results can be quite inaccurate in the case of free surface flows where a detailed sub-grid description of the interface is lacking. We address these issues by treating the large interior of the fluid region with a standard Eulerian method while (only) meshing up a thickened region near the free surface, so that both efficiency and accuracy are addressed (see Fig. 1). Moreover, with the free surface and a thickened region of the fluid near the free surface represented by a Lagrangian mesh, an exact balance between the surface tension force and the pressure force can be achieved resulting in better stability. We derive a second order accurate symmetric positive definite discretization of the Poisson equation on this hybrid particle grid structure in both two and three spatial dimensions. We then extend this Poisson equation solver to the heat equation on the cell centers, and furthermore modify the discretization to solve for the viscous forces on the MAC grid faces. Although this hybrid solver is motivated by the need for surface tension stability, it is also beneficial for simulating incompressible free surface flows without surface tension.
A common issue with particle-based solvers for incompressible flows is drifting of particles due to numerical errors, and as discussed in [19] this drifting causes variations in the particle density which subsequently cause robustness issues. In order to control particle density, various authors proposed applying a second elliptic projection [40], [36], [25]. Unlike these approaches, we calculate the particle masses from their control volumes and a constant density. In this way, the particle density is independent of the fluid density avoiding robustness issues. Although we do not require a second projection, we point out that one could combine the projection for enforcing incompressibility and the projection for enforcing a constant particle density into a single projection using the technique proposed in [38].
Section 7 addresses surface tension. For illustrative purposes, we start with a fully Lagrangian representation and we conduct the exploration using a step-by-step constructive approach. We stress that fully coupling the surface tension force to the pressure force seems to be a key for stability. In addition, imposing a constant volume condition instead of the divergence free condition seems to greatly improve accuracy especially for large time steps as in the steady state case. The final scheme is a fully implicit volume conserving scheme that can achieve a steady state solution with very large time steps. This scheme is extended to our hybrid particle grid structure for the sake of efficiency. Finally, we show that the same ideas can be applied to the standard front tracking method and to the particle level set method. Note that our implicit approach requires formulating a tightly coupled system between the implicit surface tension force and the incompressible pressure force, and we utilize the technique proposed in [42] so that this coupled system can be made symmetric positive definite and subsequently be solved efficiently.
Section snippets
Mesh generation
The entire computational domain is defined using a Marker-And-Cell (MAC) grid, and the part of the domain occupied by fluid is identified by the presence of marker particles. Particles are only needed in the thickened outer band and potentially a small buffer region, however for the sake of simplicity we utilize particles throughout the entire region occupied by fluid. We have two types of particles, surface particles and regular particles. Given an expression for the initial position and shape
Advection
Consider the advection of an arbitrary scalar field ψ driven by a given velocity field u. The ψ values are defined at the centers of the interior cells and on the particles of the Lagrangian mesh as well as the grid nodes of buffer cells. Particles that are not connected to the Lagrangian mesh do not have their ψ values readily defined, and these values can be interpolated either from the Lagrangian mesh using barycentric interpolation or from grid cells using bilinear/trilinear interpolation.
Poisson equation
Aiming at incompressible flow in Section 6, we first consider the variable coefficient Poisson equation where Ω is the interior region, and are the portions of the boundary where Dirichlet and Neumann boundary conditions are enforced, and is the outward pointing normal. We test boundary conditions similar to the incompressible flow problems of interest, i.e. Dirichlet boundary conditions everywhere, except where Ω
Heat equation
Consider the heat equation where β is the spatially constant diffusion coefficient (spatially varying viscosity is not considered). We test boundary conditions similar to what would apply to the viscosity term in the incompressible flow problems of interest, i.e. Neumann boundary conditions everywhere except on the Cartesian computational domain boundary where Dirichlet boundary conditions are applied.
We summarize the
Incompressible flow
We consider the incompressible Navier–Stokes equations as follows where ρ is the density, is the velocity, p is the pressure, ν is the kinematic viscosity, and is gravity. The entire computational domain is defined using a MAC grid, and the part of the domain occupied by fluid is discretized using a hybrid particle grid structure as described in Section 2. Velocities are stored on particles in a colocated fashion and are also stored on MAC grid faces
Surface tension
We consider incompressible Navier–Stokes equations with surface tension as follows where is the surface tension force density. For simplicity, the viscous term is ignored during the derivation in this section, but adding this force back is straightforward.
We take a step-by-step constructive approach to the discretization for illustrative purposes. We start with a fully Lagrangian representation of the fluid in Section 7.1, in light of the fact that an
Summary
We designed a new spatial discretization for simulating the incompressible Navier–Stokes equations with a free surface by hybridizing a Lagrangian mesh with the MAC grid so that the details near the interface are well resolved by particles while the majority of the interior fluid region is efficiently handled by the MAC grid. In order to accomplish this, we developed a meshing algorithm that not only correctly captures the fluid geometry but also conforms to the MAC grid in the interior. Since
Limitations and future work
There are numerous avenues for future work. In this paper, we have used marker particles throughout the entire fluid region for simplicity, and a further optimization would be to use particles only within a band near the free surface. In Section 2, the Lagrangian mesh is generated against the boundary of the computation domain, but a more efficient solution would allow MAC grid cells to touch the boundary and only mesh the region near the free surface. In Section 4, the coupled Poisson equation
Acknowledgements
Research supported in part by ONR N00014-13-1-0346, ONR N00014-09-1-0101, ONR N00014-11-1-0027, ONR N00014-11-1-0707, ARL AHPCRC W911NF-07-0027, and the Intel Science and Technology Center for Visual Computing. Computing resources were provided in part by ONR N00014-05-1-0479.
References (64)
- et al.
A monolithic mass tracking formulation for bubbles in incompressible flow
J. Comput. Phys.
(2013) - et al.
The fast construction of extension velocities in level set methods
J. Comput. Phys.
(1999) A partial differential equation approach to multidimensional extrapolation
J. Comput. Phys.
(2004)- et al.
Two-step hybrid conservative remapping for multimaterial arbitrary Lagrangian–Eulerian methods
J. Comput. Phys.
(July 2011) - et al.
FLIP: a method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions
J. Comput. Phys.
(1986) A numerical method for solving incompressible viscous flow problems
J. Comput. Phys.
(1967)- et al.
Computing a null divergence velocity field using smoothed particle hydrodynamics
J. Comput. Phys.
(2006) - et al.
An SPH projection method
J. Comput. Phys.
(1999) - et al.
Perturbations for Delaunay and weighted Delaunay 3D triangulations
Comput. Geom.
(2011) - et al.
Back and forth error compensation and correction methods for removing errors induced by uneven gradients of the level set function
J. Comput. Phys.
(2003)
A hybrid particle level set method for improved interface capturing
J. Comput. Phys.
A second-order-accurate symmetric discretization of the Poisson equation on irregular domains
J. Comput. Phys.
An incompressible multi-phase SPH method
J. Comput. Phys.
Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method
J. Comput. Phys.
Efficient simulation of surface tension-dominated flows through enhanced interface geometry interrogation
J. Comput. Phys.
A hybrid particle-mesh method for viscous, incompressible, multiphase flows
J. Comput. Phys.
Spatially adaptive techniques for level set methods and incompressible flow
Comput. Fluids
A symmetric positive definite formulation for monolithic fluid structure interaction
J. Comput. Phys.
Implicit surface tension formulation with a Lagrangian surface mesh on an Eulerian simulation grid
J. Comput. Phys.
Efficient implementation of essentially non-oscillatory shock capturing schemes
J. Comput. Phys.
An improved level set method for incompressible two-phase flows
Comput. Fluids
A front-tracking method for the computations of multiphase flow
J. Comput. Phys.
A front-tracking method for viscous, incompressible, multifluid flows
J. Comput. Phys.
Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach
J. Comput. Phys.
Fully multidimensional flux-corrected transport algorithms for fluids
J. Comput. Phys.
Efficient subdivision of finite-element datasets into consistent tetrahedra
Surface reconstruction by Voronoi filtering
Discrete Comput. Geom.
A discrete Laplace–Beltrami operator for simplicial surfaces
Discrete Comput. Geom.
A simple average nodal pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit applications
Commun. Numer. Methods Eng.
Computing Dirichlet tessellations
Comput. J.
Some improvements of the fast marching method
SIAM J. Sci. Comput.
Curve and Surface Reconstruction: Algorithms with Mathematical Analysis
Cited by (32)
Breaching the capillary time-step constraint using a coupled VOF method with implicit surface tension
2022, Journal of Computational PhysicsA point-mass particle method for the simulation of immiscible multiphase flows on an Eulerian grid
2019, Journal of Computational PhysicsCitation Excerpt :Methods that couple Lagrangian particles with an Eulerian mesh have been proposed in a variety of contexts. Lagrangian surface meshes have been used to define implicit surface tension forces in Eulerian simulations [37], and have also been used to assist in solution of the Poisson equation on hybrid Eulerian-Lagrangian meshes [38]. Particle disorganization in Lagrangian methods adversely affects the fluid pressure computed with Eulerian kernels, and the use of an Eulerian background mesh has been shown to enhance the quality of pressure computed in moving particle semi-implicit methods [39] ([39] also provides a thorough review of other efforts by researchers to use Eulerian grids, arbitrary Lagrangian-Eulerian methods, and re-mesh methods to overcome the numerical challenges associated with poorly-organized Lagrangian particles).
A semi-implicit and unconditionally stable approximation of the surface tension in two-phase fluids
2019, Journal of Computational PhysicsSharp interface approaches and deep learning techniques for multiphase flows
2019, Journal of Computational PhysicsA-SLEIPNNIR: A multiscale, anisotropic adaptive, particle level set framework for moving interfaces. Transport equation applications
2019, Journal of Computational PhysicsCitation Excerpt :Instances of interface-capturing schemes are the Volume-Of-Fluid (VOF) method [14] and improvements such as the PLIC (‘Piecewise Linear Interface Calculation’) schemes [15] or the hybrid geometric/algebraic VOF approach [16], all of them featuring excellent mass conservation properties; the Phase-Field (PF) method [17], whose diffuse interface technique offers advantages when modeling physical systems such as fracture propagation [18] or biological membranes [19]; the Immersed Boundary Method and the Level Set (LS) method [20–22], of widespread use in moving-interface problems due to its flexibility when dealing with topologically complex domains. Many other techniques build upon the previous mesh-based methods: thus, the Particle Level Set (PLS) method by Enright and co-workers [23,24], the Moment-Of-Fluid (MOF) method by Dyadechko & Shashkov [25], the hybrid particle MAC method by Zheng et al. [26], the Cellwise Conservative Unsplit VOF method by Comminal and collaborators [27], or the improved Polygonal Area Mapping (iPAM) approach by Zhang and Fogelson [28,29]. And yet, in a number of situations (thin filaments, stretching flows, boundary layers, interface merging/breaking, etc.), the spatial resolution required by mesh-based methods takes such a toll on their computational efficiency as to render them unsuitable for everyday engineering practice.
An implicit surface tension model for the analysis of droplet dynamics
2018, Journal of Computational PhysicsCitation Excerpt :In the current formulation, this factor can be as high as 200. However, this value is lower than the factor of almost 1000 reported by Zheng et al. [24]. On a recent publication [5], a two-dimensional version of the proposed model was used to predict the deformation of a sessile droplet placed on a horizontal substrate (SIGRACET 24BC, GDL side) that is then tilted at a constant rate of 0.4 deg s−1.