A new incompressibility discretization for a hybrid particle MAC grid representation with surface tension

doi:10.1016/j.jcp.2014.08.051

Journal of Computational Physics

Volume 280, 1 January 2015, Pages 96-142

https://doi.org/10.1016/j.jcp.2014.08.051 Get rights and content

Abstract

We take a particle based approach to incompressible free surface flow motivated by the fact that an explicit representation of the interface geometry and internal deformations gives precise feedback to an implicit solver for surface tension. Methods that enforce incompressibility directly on the particles are typically numerically inefficient compared to those that utilize a background grid. However, background grid discretizations suffer from inaccuracy near the free surface where they do not properly capture the interface geometry. Therefore, our incompressibility discretization utilizes a particle based projection near the interface and a background MAC grid based projection for efficiency in the vast interior of the liquid domain – as well as a novel method for coupling these two disparate projections together. We show that the overall coupled elliptic solver is second order accurate, and remains second order accurate when used in conjunction with an appropriate temporal discretization for parabolic problems. A similar second order accurate discretization is derived when the MAC grid unknowns are located on faces (as opposed to cell centers) so that Navier–Stokes viscosity can be solved for implicitly as well. Finally, we present a fully implicit approach to surface tension that is robust enough to achieve a steady state solution in a single time step. Beyond stable implicit surface tension for our novel hybrid discretization, we demonstrate preliminary results for both standard front tracking and the particle level set method.

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 $Δ t = O (Δ x^{3 / 2})$ 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 $\nabla \cdot β (\vec{x}) \nabla ϕ (\vec{x}) = f (\vec{x}), \vec{x} \in Ω$ $ϕ (\vec{x}) = g (\vec{x}), \vec{x} \in \partial Ω_{D}$ $\vec{n} \cdot \nabla ϕ (\vec{x}) = h (\vec{x}), \vec{x} \in \partial Ω_{N}$ where Ω is the interior region, $\partial Ω_{D}$ and $\partial Ω_{N}$ are the portions of the boundary where Dirichlet and Neumann boundary conditions are enforced, and $\vec{n}$ 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 $\frac{\partial ϕ (\vec{x}, t)}{\partial t} = \nabla \cdot (β \nabla ϕ (\vec{x}, t)), \vec{x} \in Ω$ $ϕ (\vec{x}, t) = g (\vec{x}, t), \vec{x} \in \partial Ω_{D}$ $\vec{n} \cdot \nabla ϕ (\vec{x}, t) = h (\vec{x}, t), \vec{x} \in \partial Ω_{N}$ 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 $\frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \nabla) \vec{u} = - \frac{1}{ρ} \nabla p + ν \nabla^{2} \vec{u} + \vec{g},$ $\nabla \cdot \vec{u} = 0,$ where ρ is the density, $\vec{u}$ is the velocity, p is the pressure, ν is the kinematic viscosity, and $\vec{g}$ 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 $\frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \nabla) \vec{u} = - \frac{1}{ρ} \nabla p + \frac{1}{ρ} \vec{f} + \vec{g},$ $\nabla \cdot \vec{u} = 0,$ where $\vec{f}$ 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)

Mridul Aanjaneya et al.
A monolithic mass tracking formulation for bubbles in incompressible flow
J. Comput. Phys.
(2013)
D. Adalsteinsson et al.
The fast construction of extension velocities in level set methods
J. Comput. Phys.
(1999)
Tariq D. Aslam
A partial differential equation approach to multidimensional extrapolation
J. Comput. Phys.
(2004)
Markus Berndt et al.
Two-step hybrid conservative remapping for multimaterial arbitrary Lagrangian–Eulerian methods
J. Comput. Phys.
(July 2011)
J.U. Brackbill et al.
FLIP: a method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions
J. Comput. Phys.
(1986)
A.J. Chorin
A numerical method for solving incompressible viscous flow problems
J. Comput. Phys.
(1967)
F. Colin et al.
Computing a null divergence velocity field using smoothed particle hydrodynamics
J. Comput. Phys.
(2006)
S. Cummins et al.
An SPH projection method
J. Comput. Phys.
(1999)
Olivier Devillers et al.
Perturbations for Delaunay and weighted Delaunay 3D triangulations
Comput. Geom.
(2011)
T. Dupont 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)

D. Enright et al.

A hybrid particle level set method for improved interface capturing

J. Comput. Phys.

(2002)

F. Gibou et al.

A second-order-accurate symmetric discretization of the Poisson equation on irregular domains

J. Comput. Phys.

(2002)

X.Y. Hu et al.

An incompressible multi-phase SPH method

J. Comput. Phys.

(2007)

E.-S. Lee et al.

Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method

J. Comput. Phys.

(2008)

P. Liovic et al.

Efficient simulation of surface tension-dominated flows through enhanced interface geometry interrogation

J. Comput. Phys.

(2010)

J. Liu et al.

A hybrid particle-mesh method for viscous, incompressible, multiphase flows

J. Comput. Phys.

(2005)

F. Losasso et al.

Spatially adaptive techniques for level set methods and incompressible flow

Comput. Fluids

(2006)

A. Robinson-Mosher et al.

A symmetric positive definite formulation for monolithic fluid structure interaction

J. Comput. Phys.

(2011)

C. Schroeder et al.

Implicit surface tension formulation with a Lagrangian surface mesh on an Eulerian simulation grid

J. Comput. Phys.

(2012)

C.-W. Shu et al.

Efficient implementation of essentially non-oscillatory shock capturing schemes

J. Comput. Phys.

(1988)

M. Sussman et al.

An improved level set method for incompressible two-phase flows

Comput. Fluids

(1998)

G. Tryggvason et al.

A front-tracking method for the computations of multiphase flow

J. Comput. Phys.

(2001)

S.O. Unverdi et al.

A front-tracking method for viscous, incompressible, multifluid flows

J. Comput. Phys.

(1992)

Rui Xu et al.

Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach

J. Comput. Phys.

(October 2009)

S.T. Zalesak

Fully multidimensional flux-corrected transport algorithms for fluids

J. Comput. Phys.

(1979)

Guy Albertelli et al.

Efficient subdivision of finite-element datasets into consistent tetrahedra

N. Amenta et al.

Surface reconstruction by Voronoi filtering

Discrete Comput. Geom.

(1999)

A.I. Bobenko et al.

A discrete Laplace–Beltrami operator for simplicial surfaces

Discrete Comput. Geom.

(2007)

J. Bonet et al.

A simple average nodal pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit applications

Commun. Numer. Methods Eng.

(1998)

A. Bowyer

Computing Dirichlet tessellations

Comput. J.

(1981)

D. Chopp

Some improvements of the fast marching method

SIAM J. Sci. Comput.

(2001)

T.K. Dey

Curve and Surface Reconstruction: Algorithms with Mathematical Analysis

(2011)

Cited by (32)

Breaching the capillary time-step constraint using a coupled VOF method with implicit surface tension
2022, Journal of Computational Physics
The capillary time-step constraint is the dominant limitation on the applicable time-step in many simulations of interfacial flows with surface tension and, consequently, governs the execution time of these simulations. We propose a fully-coupled pressure-based algorithm based on an algebraic Volume-of-Fluid (VOF) method in conjunction with an implicit linearised surface tension treatment that can breach the capillary time-step constraint. The advection of the interface is solved together with the momentum and continuity equations of the interfacial flow in a single system of linearised equations, providing an implicit coupling between pressure, velocity and the VOF colour function used to distinguish the interacting fluids. Surface tension is treated with an implicit formulation of the Continuum Surface Force (CSF) model, whereby both the interface curvature and the gradient of the colour function are treated implicitly with respect to the colour function. The presented results demonstrate that a time-step larger than the capillary time-step can be applied with this new numerical framework, as long as other relevant time-step restrictions are satisfied, including a time-step restriction associated with surface tension, density as well as viscosity.
A point-mass particle method for the simulation of immiscible multiphase flows on an Eulerian grid
2019, Journal of Computational Physics
Citation 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).
We present an Eulerian-Lagrangian approach for the modeling and simulation of immiscible multiphase flow systems. The Naiver-Stokes equations are solved on a traditional Eulerian grid while the fluid mass and phase information is discretized by Lagrangian particles. The method is novel because the particles move with a velocity that enforces consistency between the particle field density and the fluid density. The approach is advantageous in that (i) an arbitrary number of phases are easily represented, (ii) the particles remain well-distributed in space, even near merging and diverging characteristics, (iii) mass conservation is easily controlled, and (iv) the methodology is applicable to a wide range of Courant numbers. The governing equations are derived and a numerical method is presented that is applicable to incompressible flows. Performance is assessed via standard two-dimensional and three-dimensional phase transport tests as a function of both Eulerian grid resolution and Lagrangian particle resolution. Results show that the shape error converges with first-order with respect to increasing either Eulerian grid resolution or particle resolution, while mass conservation errors converge with the square root. The method is shown to successfully simulate expanding elliptical regions, stationary and oscillating droplets, a droplet in shear flow, a Rayleigh-Taylor instability, and the air blast atomization of a droplet.
A semi-implicit and unconditionally stable approximation of the surface tension in two-phase fluids
2019, Journal of Computational Physics
We focus on the energy conservation of two-phase fluids, where the change in kinetic energy is balanced by the gravitational potential and the surface energy related to surface tension. We introduce an unconditionally stable approximation in the sense that the total energy does not increase with any time step. An analysis is presented to prove the unconditional stability property of the scheme, and numerical results are given to confirming the analysis.
Sharp interface approaches and deep learning techniques for multiphase flows
2019, Journal of Computational Physics
We present a review on numerical methods for simulating multiphase and free surface flows. We focus in particular on numerical methods that seek to preserve the discontinuous nature of the solutions across the interface between phases. We provide a discussion on the Ghost-Fluid and Voronoi Interface methods, on the treatment of surface tension forces that avoid stringent time step restrictions, on adaptive grid refinement techniques for improved efficiency and on parallel computing approaches. We present the results of some simulations obtained with these treatments in two and three spatial dimensions. We also provide a discussion of Machine Learning and Deep Learning techniques in the context of multiphase flows and propose several future potential research thrusts for using deep learning to enhance the study and simulation of multiphase flows.
A-SLEIPNNIR: A multiscale, anisotropic adaptive, particle level set framework for moving interfaces. Transport equation applications
2019, Journal of Computational Physics
Citation 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.
We introduce in this paper ‘A-SLEIPNNIR’ (“Adaptive Semi-Lagrangian Ensemble Implementation of Particle level set for Newtonian and non-Newtonian Interfacial Rheology”). The proposed method combines a semi-Lagrangian approach for the transport operator, a particle level set technique for interface capturing featuring second-order accurate redistancing, and anisotropic mesh refinement for spatial resolution. The high-order, quadratic Finite Element discretization, together with an a posteriori error analysis produce a metric tensor that is passed on to an anisotropic mesh generator. We explore the peculiar features of the numerical scheme such as the impact of the marker particles, the reinitialization procedure or the Adaptive Mesh Refinement (AMR) strategy, in a series of benchmark problems for the solution of the interface transport equation with analytically supplied, divergence-free velocity fields. The results show that our approach yields an accurate and robust, yet computationally inexpensive framework for moving interfaces suitable for Scientific and Engineering applications.
An implicit surface tension model for the analysis of droplet dynamics
2018, Journal of Computational Physics
Citation 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.
A Lagrangian incompressible fluid flow model is extended by including an implicit surface tension term in order to analyze droplet dynamics. The Lagrangian framework is adopted to model the fluid and track its boundary, and the implicit surface tension term is used to introduce the appropriate forces at the domain boundary. The introduction of the tangent matrix corresponding to the surface tension force term ensures enhanced stability of the derived model. Static, dynamic and sessile droplet examples are simulated to validate the model and evaluate its performance. Numerical results are capable of reproducing the pressure distribution in droplets, and the advancing and receding contact angles evolution for droplets in varying substrates and inclined planes. The model is stable even at time steps up to 20 times larger than previously reported in literature and achieves first and second order convergence in time and space, respectively. The present implicit surface tension implementation is applicable to any model where the interface is represented by a moving boundary mesh.

View all citing articles on Scopus

View full text