Elsevier

Journal of Computational Physics

Volume 279, 15 December 2014, Pages 127-144
Journal of Computational Physics

Compressive VOF method with skewness correction to capture sharp interfaces on arbitrary meshes

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

Abstract

The accurate and efficient modelling of two-phase flows is at present mostly limited to structured, unskewed meshes, due to the additional topological and numerical complexity of arbitrary, unstructured meshes. Compressive VOF methods which discretize the interface advection with algebraic differencing schemes are computationally efficient and inherently applicable to arbitrary meshes. However, compressive VOF methods evidently suffer severely from numerical diffusion on meshes with topological skewness. In this paper we present a compressive VOF method using a state-of-the-art donor–acceptor advection scheme which includes novel modifications to substantially reduce numerical diffusion on arbitrary meshes without adding computational complexity. The new methodology accurately captures evolving interfaces on any arbitrary, non-overlapping mesh and conserves mass within the limits of the applied solver tolerance. A thorough validation of the presented methods is conducted, examining the pure advection of the interface indicator function as well as the application to evolving interfaces with surface tension. Crucially, the results on equidistant Cartesian and arbitrary tetrahedral meshes are shown to be comparable and accurate.

Introduction

The accurate representation and advection of the material interface between two fluids is essential for the predictive quality of immiscible two-phase flow simulations. However, the interface is infinitesimally thin with respect to continuum mechanics and, thus, represents a considerable challenge for finite volume and finite element frameworks. The Volume of Fluid (VOF) [1] method is among the most widely used methods to represent two incompressible, immiscible fluids. Two fundamental types of VOF methods may be distinguished: compressive methods and geometric methods. Both types of VOF methods capture the discrete volume fraction of each phase and transport it based on the underlying fluid. Compressive VOF methods discretize the partial differential equation describing the transport of the volume fraction of each phase using algebraic differencing schemes [2], [3], [4], [5]. However, the temporal and spatial discretisation requires special care in order to keep the interface sharp and without distortion. Using a geometric method, an explicit representation of the interface is advected, reconstructed from the VOF volume fraction field. The most notable state-of-the-art reconstruction methods are piecewise linear (PLIC) methods [6], [7], [8], [9], [10], [11], [12] and parabolic reconstruction methods [13]. The explicit interface is geometrically fitted to the VOF volume fraction field and advected in an Eulerian, Lagrangian or mixed Eulerian–Lagrangian fashion [14].

The major advantages of compressive VOF methods compared to geometric VOF methods are the straightforward implementation on arbitrary meshes and the computational efficiency. Geometric methods advect the interface very accurately, but the three-dimensional reconstruction of the interface requires significantly larger computational resources compared to a two-dimensional reconstruction [12], and even more so on arbitrary meshes, due to the added complexity of the geometric primitives used to reconstruct the interface. Hence, the implementation of geometric VOF methods is considerably more complex than compressive VOF methods, the required computational effort is higher and, as a results, geometric VOF methods are at present, apart from a few recently proposed exceptions [14], [15], [16], [17], exclusively available on structured meshes. Despite the greater flexibility and wider applicability of compressive VOF methods, geometric VOF methods such as PLIC methods are in current literature widely preferred over compressive VOF methods, due to their claimed higher accuracy. However, recent studies by Gopala and van Wachem [18], Park et al. [19] and Denner et al. [20] show that compressive VOF methods are generally capable of advecting sharp, evolving interfaces with similar accuracy as state-of-the-art PLIC methods.

A major challenge for compressive VOF methods is to retain the shape and sharpness of the interface. The transient advection of the interface is typically discretized using a second-order temporal discretisation scheme, such as the Crank–Nicolson scheme or the Second-Order Backward Euler scheme [4], [21], since first-order schemes are too diffusive to maintain the sharpness of the interface. Studies by Ubbink [21], Jasak [22] and Darwish [23] comprehensively demonstrate that the First-Order Backward Euler scheme as well as the First-Order Forward Euler scheme distort the shape of a circular interface advected at a constant oblique velocity on an equidistant Cartesian mesh. The Crank–Nicolson scheme, on the other hand, is able to preserve the shape of the circular interface of the same test case. In [24], Moukalled and Darwish propose a class of temporal discretisation schemes, which switch between a compressive and a high-resolution scheme based on the angle between the interface normal vector and the velocity vector. Moukalled and Darwish [24] use the Second-Order Backward Euler scheme and a compressive Euler scheme, reporting better results at medium and high Courant numbers than with other commonly used schemes, such as the Crank–Nicolson scheme.

Similarly, the careful spatial discretisation of the interface advection is essential to preserve a sharp interface. Low-order advection schemes are not suitable, as they lead to significant smearing of the interface, whereas high-order schemes result in numerical oscillations and wrinkling of the interface. Compressive VOF methods, therefore, use specifically designed spatial advection schemes based on a donor–acceptor approach. A donor (upwind) and an acceptor (downwind) cell are assigned for every cell face with respect to the underlying flow field. Based on the angle between the local interface normal vector and the normal vector of the cell face, donor–acceptor schemes blend between compressive downwind and diffusive upwind schemes. Following a donor–acceptor approach, a number of spatial discretisation schemes specifically designed for compressive VOF methods have been proposed in recent years, most notably the SURFER scheme [2], the CICSAM scheme [3], the STACS scheme [4], the HRIC scheme [5] and the HiRAC scheme [25]. The application of these schemes to arbitrary meshes is typically straightforward but is also associated with discretisation errors induced by the mesh orientation and arrangement as well as the construction of an artificial upwind node for each cell face. Zhang et al. [26] examined the influence of the upwind node extrapolation and reported a profound impact on the accuracy of the tested advection schemes. Denner and van Wachem [27] observed considerable numerical diffusion simulating the buoyancy-driven rise of a bubble on a tetrahedral mesh, suggesting the numerical diffusion is a results of mesh skewness. A more detailed examination of this issue by Denner [28] supports the assumption that mesh skewness is a major source of error with respect to the interface advection. Ubbink [21] put forward two additional possible reasons for numerical diffusion of the VOF volume fraction using donor–acceptor schemes. First, the extrapolation of the upwind node, required to determine the advection of the colour function on unstructured meshes, and, second, the implicit assumption that an interface is intersecting a cell face if both adjacent cells contain an interface. Numerical diffusion as a result of the applied differencing schemes has previously been extensively discussed in the literature [4], [21], [29], [30]. However, numerical diffusion induced by mesh skewness has, so far, not been the focus of research in the context of interface advection. This is surprising, as mesh skewness is the norm on arbitrary meshes, e.g. boundary-fitted hexahedral meshes or tetrahedral meshes, and leads to a diffusion-like error that severely affects the accuracy of spatial interpolation, as for instance the linear interpolation of face-centred values from adjacent cell-centred values required to discretize the spatial advection of the VOF volume fraction.

In this article we present a revised compressive VOF framework, based on the widely used CICSAM scheme, including modifications aiming at mitigating numerical diffusion due to mesh skewness. In particular, an implicit skewness correction is proposed which substantially decreases numerical diffusion on arbitrary meshes without notably affecting the computational complexity and mass-conservation of the VOF algorithm. It is worth mentioning that the proposed modifications are not limited to the CICSAM scheme but can also be applied to any other donor–acceptor scheme (e.g. HRIC or STACS). The comprehensive and thorough validation demonstrates the capabilities and accuracy of the presented methodology and the results on Cartesian and tetrahedral meshes are shown to be comparable. Moreover, we show that our framework is mass-conserving within the limits of the applied solver tolerance. The presented methodology including the proposed modifications are of significant interest to the two-phase flow modelling community as it improves the predictive quality of interfacial flow simulations considerably and, thus, makes it feasible for applications with complex geometries which cannot be represented by general structured meshes.

The article is structured as follows. In Section 2 the governing equations are outlined and Section 3 briefly explains the applied numerical framework. The discretisation of the applied compressive VOF method is discussed in Section 4. Subsequently, Section 5 is concerned with numerical diffusion caused by mesh skewness and introduces the proposed modifications to mitigate numerical diffusion. The presented compressive VOF methodology including the proposed modifications are comprehensively tested and validated in Section 6. The findings are summarised and the article is concluded in Section 7.

Section snippets

Governing equations

The Volume of Fluid (VOF) method [1] is used to determine the position and behaviour of the interface between two incompressible, isothermal and immiscible fluids. The VOF method applies a volume fraction γ, also called colour function, to every mesh cell, representing the local volume fraction asγ(x,t)={0fluid A1fluid B. Thus, every mesh cell where the colour function value is between 0 and 1 contains an interface. As the interface is transported with the underlying flow field, the colour

Numerical framework

The numerical framework follows a coupled implicit approach for the fluid flow with a segregated interface advection, as presented by Denner and van Wachem [27]. A linear equation system describing the flow is constructed, containing the discretized momentum equations as well as a fourth equation, which provides an additional relationship between pressure and velocity to close the equation system. The fourth equation of the flow equation system is constructed based on the continuity equation,

Compressive VOF method

In the following sections, the discretisation of the interface advection equation is presented and discussed. The numerical schemes presented in this section represent only one specific choice, deemed to be best suited for the discretisation of the interface advection on arbitrary meshes. Other schemes may be used in a similar fashion without affecting the applicability of the proposed modifications presented in Section 5.

Skewness correction

Skewness is a commonly encountered complication for non-Cartesian, arbitrary meshes, as for instance depicted in Fig. 2. As a result of skewness, the geometric face centre does not coincide with the interpolation point at the face, as for instance the point where the vector connecting the adjacent cell centres intersects the face, denoted in Fig. 2 as f. Thus, the linearly interpolated face value becomes inaccurate and the accuracy of the interpolation reduces formally to first order. As

Results

The compressive VOF methodology described in Section 4, as well as the remedies for mesh skewness proposed in Section 5 are validated and compared using three representative test cases. The pure advection of a fluid particle in a constant velocity field, see Section 6.1, and in a shear flow, see Section 6.2, is simulated. These test cases allow the assessment of the proposed methods regarding accuracy and mass conservation without the influence of velocity gradients, pressure gradients or

Conclusion

In this article, the discretisation and implementation of a compressive VOF method for three-dimensional arbitrary meshes has been presented. The presented compressive VOF methodology is applicable to arbitrary, unstructured meshes and conserves the mass of each phase within the limits of the applied solver tolerance. However, as reported in previous studies [4], [21], [27], [28], the predictive quality of compressive VOF methods on unstructured meshes is severely challenged by numerical

Acknowledgement

The authors are grateful to the Engineering and Physical Sciences Research Council (EPSRC) for their financial support (grant EP/K008595/1).

References (48)

  • Z. Wang et al.

    A new volume-of-fluid method with a constructed distance function on general structured grids

    J. Comput. Phys.

    (2012)
  • V.R. Gopala et al.

    Volume of fluid methods for immiscible-fluid and free-surface flows

    Chem. Eng. J.

    (2008)
  • F. Denner et al.

    Comparative study of mass-conserving interface capturing frameworks for two-phase flows with surface tension

    Int. J. Multiph. Flow

    (2014)
  • Y.-Y. Tsui et al.

    Flux-blending schemes for interface capture in two-fluid flows

    Int. J. Heat Mass Transf.

    (2009)
  • Y. Morinishi et al.

    Fully conservative higher order finite difference schemes for incompressible flow

    J. Comput. Phys.

    (1998)
  • J.U. Brackbill et al.

    Continuum method for modeling surface tension

    J. Comput. Phys.

    (1992)
  • B.P. Leonard

    The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection

    Comput. Methods Appl. Mech. Eng.

    (1991)
  • I. Demirdžić et al.

    Numerical method for coupled fluid flow, heat transfer and stress analysis using unstructured moving meshes with cells of arbitrary topology

    Comput. Methods Appl. Mech. Eng.

    (1995)
  • M. Pivello et al.

    A fully adaptive front tracking method for the simulation of two phase flows

    Int. J. Multiph. Flow

    (2014)
  • J. Hua et al.

    Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method

    J. Comput. Phys.

    (2008)
  • F. Denner et al.

    On the convolution of fluid properties and surface force for interface capturing methods

    Int. J. Multiph. Flow

    (2013)
  • P. Lubin et al.

    Three-dimensional large eddy simulation of air entrainment under plunging breaking waves

    Coast. Eng.

    (2006)
  • M. Darwish et al.

    Convective schemes for capturing interfaces of free-surface flows on unstructured grids

    Numer. Heat Transf., Part B, Fundam.

    (2006)
  • S. Muzaferija et al.

    A two-fluid Navier–Stokes solver to simulate water entry

  • Cited by (68)

    • Modeling of two-phase flows at low Capillary number with VoF method

      2023, Computers and Fluids
      Citation Excerpt :

      To avoid all the problems connected to the numerical setup for the interface, geometrical methods reconstruct the interface by generating planes with the orientation extrapolated by the indicator function, like the PLIC (piecewise linear interface calculation) [3,19,20] and the SLIC (simple line interface calculation) [21,22] methods. The great advantage of this geometrical treatment of the interface region is an increase of accuracy [20], but they are more complex to manage and less versatile with unstructured grid [23]. In the next section, the methodology is presented, describing all its features and the CICSAM scheme definition.

    View all citing articles on Scopus
    View full text