Compressive VOF method with skewness correction to capture sharp interfaces on arbitrary meshes
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 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 . 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)
- et al.
Volume of fluid (VOF) method for the dynamics of free boundaries
J. Comput. Phys.
(1981) - et al.
Modelling merging and fragmentation in multiphase flows with SURFER
J. Comput. Phys.
(1994) - et al.
A method for capturing sharp fluid interfaces on arbitrary meshes
J. Comput. Phys.
(1999) - et al.
Reconstructing volume tracking
J. Comput. Phys.
(1998) - et al.
Second-order accurate volume-of-fluid algorithms for tracking material interfaces
J. Comput. Phys.
(2004) - et al.
Interface reconstruction with least-squares fit and split advection in three-dimensional Cartesian geometry
J. Comput. Phys.
(2007) - et al.
Analytical and geometrical tools for 3D volume of fluid methods in general grids
J. Comput. Phys.
(2008) - et al.
Simulations of multidimensional interfacial flows by an improved volume-of-fluid method
Int. J. Heat Mass Transf.
(2013) - et al.
PROST: a parabolic reconstruction of surface tension for the volume-of-fluid method
J. Comput. Phys.
(2002) - et al.
A PLIC–VOF method suited for adaptive moving grids
J. Comput. Phys.
(2011)
A new volume-of-fluid method with a constructed distance function on general structured grids
J. Comput. Phys.
Volume of fluid methods for immiscible-fluid and free-surface flows
Chem. Eng. J.
Comparative study of mass-conserving interface capturing frameworks for two-phase flows with surface tension
Int. J. Multiph. Flow
Flux-blending schemes for interface capture in two-fluid flows
Int. J. Heat Mass Transf.
Fully conservative higher order finite difference schemes for incompressible flow
J. Comput. Phys.
Continuum method for modeling surface tension
J. Comput. Phys.
The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection
Comput. Methods Appl. Mech. Eng.
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.
A fully adaptive front tracking method for the simulation of two phase flows
Int. J. Multiph. Flow
Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method
J. Comput. Phys.
On the convolution of fluid properties and surface force for interface capturing methods
Int. J. Multiph. Flow
Three-dimensional large eddy simulation of air entrainment under plunging breaking waves
Coast. Eng.
Convective schemes for capturing interfaces of free-surface flows on unstructured grids
Numer. Heat Transf., Part B, Fundam.
A two-fluid Navier–Stokes solver to simulate water entry
Cited by (68)
A novel steepness-adjustable harmonic volume-of-fluid method for interface capturing
2024, Journal of Computational PhysicsAn accurate and efficient multiphase solver based on THINC scheme and adaptive mesh refinement
2023, International Journal of Multiphase FlowModeling of two-phase flows at low Capillary number with VoF method
2023, Computers and FluidsCitation 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.
Effects of cell quality in grid boundary layer on the simulated flow around a square cylinder
2022, Computers and FluidsModeling interfacial mass transfer of highly non-ideal mixtures using an algebraic VOF method
2022, Chemical Engineering ScienceA grouting simulation method for quick-setting slurry in karst conduit: The sequential flow and solidification method
2022, Journal of Rock Mechanics and Geotechnical Engineering