A novel consistent and well-balanced algorithm for simulations of multiphase flows on unstructured grids

Journal of Computational Physics 350 (2017) 207–236

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A novel consistent and well-balanced algorithm for simulations of multiphase flows on unstructured grids Jitendra Kumar Patel, Ganesh Natarajan ∗ Department of Mechanical Engineering, Indian Institute of Technology Guwahati, Assam 781039, India

a r t i c l e

i n f o

Article history: Received 21 January 2017 Received in revised form 11 June 2017 Accepted 23 August 2017 Available online xxxx Keywords: Volume-of-fluid Consistent transport Balanced force algorithm High density ratios Hybrid staggered/non-staggered Multiphase flows

a b s t r a c t We discuss the development and assessment of a robust numerical algorithm for simulating multiphase flows with complex interfaces and high density ratios on arbitrary polygonal meshes. The algorithm combines the volume-of-fluid method with an incremental projection approach for incompressible multiphase flows in a novel hybrid staggered/non-staggered framework. The key principles that characterise the algorithm are the consistent treatment of discrete mass and momentum transport and the similar discretisation of force terms appearing in the momentum equation. The former is achieved by invoking identical schemes for convective transport of volume fraction and momentum in the respective discrete equations while the latter is realised by representing the gravity and surface tension terms as gradients of suitable scalars which are then discretised in identical fashion resulting in a balanced formulation. The hybrid staggered/non-staggered framework employed herein solves for the scalar normal momentum at the cell faces, while the volume fraction is computed at the cell centroids. This is shown to naturally lead to similar terms for pressure and its correction in the momentum and pressure correction equations respectively, which are again treated discretely in a similar manner. We show that spurious currents that corrupt the solution may arise both from an unbalanced formulation where forces (gravity and surface tension) are discretised in dissimilar manner and from an inconsistent approach where different schemes are used to convect the mass and momentum, with the latter prominent in flows which are convection-dominant with high density ratios. Interestingly, the inconsistent approach is shown to perform as well as the consistent approach even for high density ratio flows in some cases while it exhibits anomalous behaviour for other scenarios, even at low density ratios. Using a plethora of test problems of increasing complexity, we conclusively demonstrate that the consistent transport and balanced force treatment results in a numerically stable solution procedure and physically consistent results. The algorithm proposed in this study qualifies as a robust approach to simulate multiphase flows with high density ratios on unstructured meshes and may be realised in existing flow solvers with relative ease. © 2017 Elsevier Inc. All rights reserved.

1. Introduction The interaction of immiscible fluids with distinct interface(s) can be seen in several engineering and environmental phenomena. Understanding these multi-component flows and their underlying physics necessitates the use of high-resolution

*

Corresponding author. E-mail addresses: [email protected] (J.K. Patel), [email protected] (G. Natarajan).

http://dx.doi.org/10.1016/j.jcp.2017.08.047 0021-9991/© 2017 Elsevier Inc. All rights reserved.

208

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

numerical algorithms for their simulations. Unlike incompressible approaches for single-phase flows, studies in multiphase flows are fraught with challenges that include resolving the interface(s) sharply and simulating realistic and high density and viscosity contrasts successfully. In addition, the presence of interfacial forces due to surface tension and volumetric forces such as those due to gravity and electromagnetic effects need to be suitably accounted for in the numerical discretisation to ensure a stable solution procedure. Multiphase simulations have been carried out using several different approaches over the past two decades. While a comprehensive and exhaustive review of the various methods may be found in [1], we briefly recall the salient contributions in this field of research. Restricting to a single fluid formalism, it is easy to see that the widely used approaches belong to the class of volume-of-fluid (VoF) or the level-set methods with both approaches having their distinct strengths and weaknesses. The VoF approach, which originated from the work of Hirt and Nichols [2] while being discretely conservative leads to discontinuous property variations whereas the level-set approach pioneered by Osher and Sethian [3] results in a smoother transition of fluid properties at the cost of sacrificing conservation. There have been attempts at developing variants such as the mass-conserving level set method (MCLS) [4] and also at hybridising these approaches such as the CLSVOF method [5] and the VOSET method [6], which have met with considerable success. Simulations of large density multiphase flows have also found significant attention in the framework of lattice Boltzmann methods [7] as well as diffuse interface approaches [8]. An important issue that plagues multiphase flow simulations, particularly those involving interfacial tension, is the generation of spurious currents. Spurious currents in surface-tension dominated flows arise due to the errors in calculation of curvature and the discrete imbalance between surface tension and pressure forces [9]. Curvature errors can be minimised using improved approaches for calculating the interface curvature as detailed in [10–12] but can never be completely eliminated and represent a constant source of non-zero spurious currents. Abadie et al. [13] have shown that spurious currents are related to spurious vorticity production and have investigated the role of advection schemes in their evolution in conjunction with different curvature calculation approaches. It is also discussed in [14] that coupling between surface tension force, advection scheme and momentum transport is crucial. A recent study by Pan et al. [15] have shown that spurious currents in low Capillary number flows can be significantly reduced by employing a moving reference frame for simulations. However, one can totally eliminate errors due to force imbalances by carefully constructing a discretisation that enforces balance between competing forces, leading to a class of approaches referred to as “well-balanced” (or “balanced force”) algorithms. The principle behind these class of methods is to use identical discretisation of the pressure gradient and surface tension term that appear in the momentum equation and with a reasonably accurate interface curvature calculation leads to acceptably low spurious currents. The earliest and significant contributions to well-balanced algorithms include the works of Renardy and Renardy [16] and Francois et al. [17]. While the former proposed the parabolic reconstruction of surface tension (PROST) for accurate and balanced implementation of interfacial tension, the latter discusses the need for and development of a balanced algorithm for multiphase flows on uniform structured meshes. Their algorithm however involved interpolations of centroidal quantities to faces and vice-versa and did not account for body forces. Montazeri and Ward [18] devised an improved pressure-velocity coupling based on Newton’s divided-differences for multiphase flows that is effective on non-uniform meshes and considers interfacial and body forces. Denner and van Wachem [19] have proposed a pressure-velocity coupling for multiphase flows on unstructured meshes with modifications to the original momentum interpolation method for surface tension and gravity forces as well as grid non-orthogonality. Ghods and Herrmann [20] adopt a different approach to achieve the force balance by employing a least-squares based reconstruction for centroidal gradients, although despite its generic nature this approach has not been extensively investigated on unstructured grids for multiphase flows with large body forces. It must be remarked that while many researchers discuss and test balanced force algorithms for interfacial tension dominated flows, they are equally important in presence of other forces as well, such as body forces like gravity [21,22]. A related pertinent challenge that has received relatively lesser attention is the ability of the flow solver to handle large density and viscosity ratios. This is critical because the density and viscosity ratios are O(100)–O(1000) for engineering applications and the key idea for a robust and stable numerical algorithm at high density ratios is to ensure a consistent transport of mass and momentum and the discrete level. This was first proposed by Rudman within a VoF framework on staggered grids [23] which is also adopted to level-set method by Desjardins and Moureau [24] and later extended to unstructured collocated meshes by Bussmann [25]. Raessi and Pitsch [26] proposed a density-based momentum flux correction in a level-set framework for tightly coupling mass and momentum, where the density fluxes computed from the level-set are used to define the momentum fluxes. More recently, Ghods and Herrmann [20] proposed a consistent rescaled mass-momentum transport (CRMT) algorithm in a level-set framework for multiphase flows. The fundamental philosophy was to convect mass and momentum using the same upwind scheme, then solve the level set equation using a higher-order accurate scheme and recompute the density using the level set information. This led to conservation errors in both mass and momentum, but in a consistent manner, which allowed a stable simulation of large density ratio flows. The focus of the present work is to develop and investigate a fully balanced and consistent transport algorithm for multiphase flows on arbitrary polygonal meshes. The work employs a novel staggered/non-staggered framework for solving the governing equations and will be shown to naturally lead to a discretisation that motivates a balanced force approach which accounts for interfacial as well as body forces. Following Ghods and Herrmann [20], we implement a consistent mass-momentum transport in the VoF framework with a second-order accurate bounded scheme that is shown to work

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

209

for arbitrarily large density ratios even in flows with low Capillary and Froude numbers. We perform extensive numerical experiments in surface-tension dominated and gravity-dominated regimes to highlight the importance and utility of our implementation on structured and unstructured grids. The remainder of the paper is organised as follows. The governing equations for the multiphase flow and the numerical methodology are presented in Sections 2 and 3 respectively. Section 4 discusses the implementation of balanced force approach in a novel hybrid staggered/non-staggered framework. In Section 5, we describe the implementation of consistent mass and momentum transport at the discrete level. We investigate the potential of the consistent and well balanced algorithm through a series of numerical experiments in Section 6. In Section 7, we briefly discuss the implications of consistent and inconsistent formulations for multiphase flows in view of the results from numerical experiments. Subsequently, in Section 8, we present a comparison of our algorithm with other alternative multiphase approaches to highlight the novelty of the proposed algorithm. Finally, Section 9 summarises the findings and concludes the present study. 2. Governing equations The incompressible laminar fluid flow in two space dimensions is governed by the Navier–Stokes equations, which in conservative form can be written as,

∇ ·u=0 ∂(ρ u) + ∇ · (ρ uu) = −∇ p´ + ∇ · (μ(∇ u + ∇ uT )) + ρ g + Fst + Fb ∂t

(1) (2)

where u represents the velocity vector, g = [0, − g ] defines the gravitational acceleration vector with g being its magnitude, μ is the dynamic viscosity and ρ represents density. The terms Fst and Fb represent the forces due to surface tension and body forces (other than due to gravity). Following [18], we rewrite this equation by defining the piezometric pressure as, p = p´ − ρ (g · x), where x is the position vector. The use of piezometric pressure results in gravity being treated as an interfacial force rather than a volumetric one. Consequently, we see that both pressure and gravity forces appear as gradients of scalar quantities,1 which allows for a balanced force formulation as in [18],

∂(ρ u) + ∇ · (ρ uu) = −∇ p + ∇ · (μ(∇ u + ∇ uT )) − g · x∇ ρ + Fst + Fb ∂t

(3)

In our studies, we consider two immiscible fluids with different densities and viscosities (in general) and nondimensionalise the governing Navier–Stokes equations using the properties of less denser fluid (defined with subscripts L) as reference quantities along with suitable length and velocity scales. The resulting momentum equation then reads (with all the quantities now being dimensionless),

∇ρ ∂(ρ u) 1 + ∇ · (ρ uu) = −∇ p + ∇ · (μ(∇ u + ∇ uT )) − e · x + Fˇ st + Fˇ b ∂t Re Fr

(4)

ρL U D U2 and F r = are the Reynolds number and Froude number respectively with D and U being suitably μL gD

where Re =

defined length and velocity scales, depending on the flow problem. Furthermore, e = [0, 1] represents the unit vector in the vertical direction (parallel to gravity). The term, Fˇ st , is the dimensionless force due to surface tension, which is implemented using the continuum surface force (CSF) model [27] as,

Fˇ st =

κ nI δ

where W e =

We

,

(5)

ρL U 2 D is the Weber number with σ denoting the surface tension coefficient. Here κ , n I and δ represent the σ

interface curvature, normal and the Dirac delta function, respectively. In the VoF framework, adopted in the present study, the volume fraction φ is convected with the local flow velocity and can be used to compute the interface normal and curvature. The equation of the volume fraction reads,

Dφ Dt

=

∂φ ∂φ + u · ∇φ = + ∇ · (uφ) = 0 ∂t ∂t

(6)

where the continuity constraint has been used to rewrite the convective term in a conservative form. The curvature appearing in the computation of surface tension force is related to the interface normal as,

κ = −∇ · n I 1

Every vector may be written as the gradient of a scalar and one can easily show ∇(ρ (g · x)) = (g · x)∇ ρ + ρ ∇(g · x) = (g · x)∇ ρ + ρ g.

(7)

210

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

Fig. 1. The control volume  = m ∪ n for the normal momentum that combines the two adjacent cells m and n shared by the face f .

Following [28] we use the smoothed volume fraction to calculate the curvature and its cell-averaged value is computed as,

κ =−

1 

c



f

∇ φ˜ ˜ |∇ φ|

 · s f

(8)

where c and s f are the cell volume and surface area vector. In order to reduce errors in curvature calculations, we employ smoothed volume fractions φ˜ obtained using Laplacian smoothing [28] given by,



φ f | s f |

f

φ˜ = 

(9)

| s f |

f

where the summation is over the faces of a cell. In this work we typically use 3–4 smoothing iterations for curvature calculations while the non-smoothed volume fractions are employed everywhere else when required. Expressing the interface normal in terms of the volume fraction, Eq. (5) can be simplified as,

Fˇ st =

κ ∇φ We

(10)

.

The properties of heavy and light fluids (with subscripts H and L) in a single fluid formalism are then calculated using a volume fraction weighted averaging as,

ρh ρl μh μ = φ + (1 − φ) μl ρ = φ + (1 − φ)

(11) (12)

3. Numerical methodology The governing equations described in the earlier section are discretised using the finite volume approach employing a novel hybrid staggered/non-staggered framework. This framework provides for a “natural” extension of the staggered philosophy to unstructured meshes and proves beneficial for multiphase simulations as will be demonstrated later in the manuscript. The key component of the framework is the solution of a single momentum equation for the normal momentum at the cell faces, as opposed to individual velocity components at cell centres. The scalars (pressure and volume fraction in this study) are solved for, similar to a collocated framework, at the cell centres. The normal momentum equation is obtained by projecting Eq. (4) in the direction normal to the cell face. Integrating the normal momentum equation over the control volume  = m ∪ n (See Fig. 1) and using the Gauss divergence theorem and single-point quadrature gives the following discrete form of the Eq. (4).

δ(ρ f U f ) δt

+

  δ p n 1 

(ρe ue U e ) se · n f = −   + (μe ∇(ue + uT e ) · ne ) se · n f δn f Re e ∈ e ∈  

 



CONV

VISC

κ f δφ n+1 δ ρ n+1 − (e · x) f  +   + (Fˇ bf · n f )n+1  Fr δn f W e δn f 1

(13)

In this equation ne and se represent the outward unit normal and surface area of the edge ‘e’ respectively while the superscript denotes the solution at an intermediate time level, which represents the predictor phase. The derivaδρ δφ δp tives δn | f , δn | f and δn | f indicate the discrete form of the normal gradients of pressure, density and volume fraction at the face f and their specific implementation are discussed later in the manuscript. It needs to be noted that while the pressure gradients are at time level n, the density and volume fraction gradients are at time level (n + 1) because

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

211

the volume fraction equation is solved before the normal momentum equation. The terms CONV and VISC collectively represent the convective and viscous fluxes, with the face values of density and viscosity obtained using volumetric and harmonic averaging respectively. The temporal discretisation is effected using a three-point backward differencing scheme which reads,

δ(ρ f U f ) δt

=

ρ nf +1 (3U f − 4U nf + U nf −1 ) 2 t

(14)

Employing an incremental fractional step approach similar to that in [29], but modified to account for variable density flows, results in a variable coefficient Poisson equation for pressure correction (P = pn+1 − pn ),



U f s f =

2 t 

f

3

f

1

ρ

n +1 f

δP   s f δn f

(15)

where the summation is over faces of the control volume. It is important to remark that the control volumes for normal momentum and pressure correction equations are clearly different in the present algorithm. While the control volume for pressure are the cells themselves, that for normal momentum is the union of two cells sharing the face. However, the calculation of convective and viscous fluxes are carried out in the same manner as is done on a collocated mesh (hence also the non-staggered aspect) with the latter obtained using central differencing and the former with a suitable upwind-biased scheme. The centroidal velocities are then necessary for the these flux computations and are obtained from the face normal velocities using a special reconstruction procedure [30]. A slightly detailed derivation of the above equations in the hybrid staggered/non-staggered framework is presented in Appendix A. The normal momentum and pressure are then advanced in the corrector phase as,

p n +1 = P + p n

2 t δP  (ρ U f )n+1 = (ρ U f ) −  3 δn f

(16) (17)

It must be noted that we choose ρ = ρ n+1 since the volume fraction equation is first solved before the predictor– corrector phases in the present algorithm and density follows from Eq. (11). The volume fractions are obtained by solving the discrete volume fraction equation that reads,

c

3φcn+1 − 4φcn + φcn−1 2 t

+



φ nf +1 U nf s f = 0

(18)

f

where the same temporal discretisation as for the momentum equation is employed. The convective scheme used for volume fraction advection has been discussed in a separate work [31], where it was shown to that conserve the mass quite accurately. The algebraic volume-of-fluid approach adopted herein is free of any interface reconstruction and a brief exposition of this interface capturing approach is presented in Appendix B. The non-linear scalar equation for normal momentum is solved using a Newton–Krylov solver (through the PetSc libraries [32]), while the system of linear algebraic equations resulting from the Poisson equation and volume fraction equation are solved using SAAMG-preconditioned GMRES and ILU-preconditioned GMRES solvers respectively, employing the LiS libraries [33]. The error tolerances for all linear iterative solvers in this work were set to a low value of 10−12 . We remark that the flow solver is spatially and temporally second order accurate on arbitrary polygonal meshes but defer from discussing further details because it does not constitute a focal point of this study in itself and also because the basic framework and its assessment are part of an independent publication. The overall solution methodology of the proposed hybrid staggered/non-staggered framework for multiphase flows is summarised below.

Algorithm 1: Hybrid staggered/non-staggered solver for Navier–Stokes equations. 1. 2. 3. 4. 5. 6. 7. 8.

Initialise the velocities (u , v ), pressure ( p ) and volume fraction (φ). Solve the advection equation (Eq. (18)) for volume fraction. Update the density (ρ n+1 ) and viscosity (μn+1 ) using Eqs. (11) and (12). Solve the normal momentum equation (Eq. (13)) for the face normal momentum (ρ n+1 U f ). Solve the Poisson equation (Eq. (15)) for pressure correction (P). Update pressure ( pn+1 ) and the face normal momentum (ρ U f )n+1 following Eqs. (16) and (17). Reconstruct the centroidal velocities (un+1 , v n+1 ) from the face normal velocity. Go to step 2 and repeat the process for a finite number of outer-iterations (for an unsteady problem) or until convergence (for a steady-state problem).

212

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

Fig. 2. The projections of the centroids on normal to the face f .

4. Balanced force algorithm As discussed in the introduction, one of the critical aspects of numerical discretisation for multiphase flows is the need to ensure a balance of forces in the momentum equation. While collocated approaches in [17–19] adopt different methodologies that require modifications in momentum interpolation to achieve the balance, the present staggered/non-staggered framework allows for a relatively straightforward and arguably “natural” approach for well-balanced discretisation. A wellbalanced algorithm is one that achieves a discrete balance between forces leading to accurate predictions of stationary solutions on any mesh topology. One can easily see that the terms corresponding to pressure, surface tension and body forces in the normal momentum equation largely involve gradients of scalar quantities viz. p, φ and ρ respectively. The key to a well-balanced formulation is to discretise all these gradients in an identical manner. Unlike the collocated framework, the forces in the normal momentum equation involve only the normal derivatives (see Eq. (13)) which may be discretised accounting for non-orthogonal effects as,

∇ αn · rn − ∇ αm · rm δ α  αn − αm +  = δn f   r

  r

orthogonal part

(19)

non-orthogonal part

where α is any scalar quantity (= ρ , p , φ, P) and r denotes the distance between the points m and n . The quantities rm , rn are the vectors along the line joining the points m and m , n and n respectively (See Fig. 2) and contribute to the non-orthogonal correction. It must also emphasised that normal derivatives of pressure in the momentum equation, Eq. (13) and momentum update, Eq. (17) must necessarily have the same discretisation as also the pressure correction in the Poisson equation, Eq. (15). The gradients appearing in the non-orthogonal contribution are estimated using a standard Green–Gauss approach [34]. A different and possibly naive approach to determine the normal gradients appearing in the force terms is to use an alternate discretisation defined by,

δ α  δα δα nx, f + n y, f  = δn f δx δy

(20)

where the quantities δδαx and δδαy are obtained using a linear interpolation of these gradients at centroids of cells (m and n) sharing the face ‘ f ’ (See Fig. 2) and nx, f and n y , f are the components of the unit normal to the face. The use of Eq. (19) for pressure gradient and Eq. (20) for density and/or volume fraction gradients (or vice-versa) would however result in an unbalanced formulation. We shall discuss its implications for surface tension, electromagnetic and gravity dominated flows in Section 6. 5. Consistent discretisation of governing equations The transport of mass and momentum in a consistent fashion at the discrete level is necessary to obtain robust and stable solutions for high density and high shear flows. In the present work, we adopt the philosophy in [20] within a VoF framework to achieve discrete consistency for the transport equations on arbitrary polygonal meshes. The inherent and fundamental principle behind a consistent approach is to ensure that the transporting velocity and transported quantities are discretely similar in both the mass and momentum equations. This means that the convective fluxes in the VoF equation, Eq. (18), as well as the normal momentum equation, Eq. (13), must be evaluated necessarily using the same convective scheme. In this study, we employ a high-resolution upwind-biased scheme referred to as CUI (Cubic Upwind Interpolation) [31, 35] for the discretisation of the convective fluxes in the momentum and volume fraction equations. One can obtain the

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

213

Table 1 Spurious currents and pressure jump errors for inviscid static droplet with exact curvature on structured and unstructured mesh. The timestep is fixed to t = 0.001 with W e = 1. Time

Unbalanced

t 50 t

Balanced

t 50 t

ρh ρl

Structured

Unstructured

L ∞ (u)

E ( ptotal )

L ∞ (u)

E ( ptotal )

10 103 10 103

5.41 × 10−6 9.73 × 10−6 2.11 × 10−3 7.11 × 10−3

1.47 × 10−1 2.43 × 10−1 4.39 × 10−1 4.96 × 10−1

2.97 × 10−5 4.47 × 10−5 9.84 × 10−3 1.01 × 10−2

7.56 × 10−1 8.52 × 10−1 8.41 × 10−1 8.95 × 10−1

10 103 10 103

1.09 × 10−14 1.89 × 10−14 1.71 × 10−13 2.01 × 10−13

1.15 × 10−3 1.20 × 10−3 2.10 × 10−3 2.16 × 10−3

3.95 × 10−12 6.21 × 10−12 3.32 × 10−11 6.32 × 10−11

8.52 × 10−3 8.73 × 10−3 8.16 × 10−3 9.24 × 10−3

normalised face values of quantities in the convective fluxes using this scheme as,

⎧ ˜C 3α ⎪ ⎪ ⎪ ⎨ 5 α˜ + α˜f = 6 C ⎪1 ⎪ ⎪ ⎩ α˜ C where

˜C ≤ 0<α 1 3

2 13 ˜ C ≤ 45

2 < 13 4 ˜C < 5

α α ≤1

(21)

elsewhere

α is any transport quantity (= u , v , ρ , φ) and α˜ =

α − αU , with the subscripts C , D and U referring to the upα D − αU

wind, downwind and far upwind cells respectively for the face f . The CUI scheme is nominally second-order accurate on unstructured meshes and ensures boundedness by construction [31]. It is possible to use other high-order accurate schemes as long as the same scheme is employed for convection of volume fraction and momentum. The use of different schemes for volume fraction and momentum transport leads to an inconsistent approach. We show that an inconsistent transport, realised herein by using a first-order upwind scheme in place of CUI for momentum transport could be detrimental for multiphase flows under some conditions as discussed in Section 6. It must be also be remarked that any choice of dissimilar schemes for mass and momentum transport would also lead to inconsistency, irrespective of their orders of accuracy as demonstrated later in the manuscript. 6. Numerical tests

We investigate the role and need of consistent and fully-balanced formulations, as proposed in Sections 4 and 5 for multiphase simulations through a series of numerical experiments in this section. A wide range of test problems of varying levels of complexity are tested on both structured and unstructured meshes to analyse the implications of unbalanced formulation as well as inconsistent transport. These studies are then employed to arrive at salient conclusions on the performance of the proposed algorithm. We emphasise that the all dimensionless numbers in the studies are defined based on the properties of the less denser fluid (denoted by subscript L) which is also considered to be the less viscous one. 6.1. Static droplet We first consider the case of a static droplet in an inviscid fluid in the absence of gravity. The droplet of diameter D = 0.5 is placed at the centre of a 1 × 1  square cavity. It must be remarked that the velocity scale used in defining the

dimensionless numbers is chosen here as

σ ρL D , while the length scale is the diameter of the droplet. Studies are carried out

on a 50 × 50 Cartesian grid as well as a triangulated grid with 2030 elements at two different density ratios. All simulations are carried out at W e = 1 with a timestep t = 0.001. The balanced and unbalanced formulations are compared by considering the following two quantities,

L ∞ (u) = max|u| E ( ptotal ) =

| pn − Wκ e | κ

We

While the first is a measure of spurious currents, the second expresses the deviation of the numerically computed pressure jump ( pn = p in − p out , difference between the average computed pressures inside and outside of the droplet) from the exact pressure jump as given by the Young–Laplace law [17]. Table 1 shows the spurious currents and pressure error jumps on the two grids when exact curvature is specified. It is evident that the balanced algorithm maintains near machine-zero values of spurious currents on either mesh even at large times, independent of the density ratio. On the contrary, the unbalanced formulation produces comparatively larger

214

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

ρ

Fig. 3. The pressure errors E ( ptotal ) in static inviscid droplet showing (a) spatial and (b) temporal accuracies for ρh = 1000 and W e = 1 using exact l curvature on unstructured grid. The dashed line represents the second order accuracy.

Fig. 4. The evolution of TKE with time for inviscid static droplet using Ubbink curvature model in (a) unbalanced and (b) balanced algorithms. The simulaρ tions are performed for density ratio ρh = 10 and W e = 1 with timestep t = 0.001. l

currents which quickly increase with time. The effect of the balanced algorithm can also be observed in the errors in pressure jump where the errors are lower by two orders of magnitude as compared to the unbalanced formulation. It must be remarked that the errors are larger on unstructured grids as compared to the structured grids for both formulations and can be attributed to the relative coarseness of the unstructured grids employed. However, the errors in pressure jump on the structured meshes using the balanced formulation in the present study are lower than those reported in a similar study [17], where a similar grid resolution was employed. This may be possibly attributed to differences in the implementation of the balanced algorithm between the present study and those in [17]. We repeat this test case for a density ratio equal to 1000 to study the spatio-temporal accuracy on unstructured meshes. The spatial accuracy is investigated by considering four different progressively refined unstructured meshes (starting from 2030 triangular elements) with a constant timestep t = 0.001. In order to obtain the temporal accuracy, the study was carried out using four different timesteps ( t = 0.005, 0.0025, 0.00125, 0.000625) on the coarsest grid. It can be seen in Figs. 3(a) and (b) that the proposed algorithm shows nominal second-order accuracy both in space and time, making it an accurate flow solver for multiphase applications. We further consider the performance of the two formulations when the interface curvature is calculated. The evolution of the total kinetic energy, which must be ideally zero, is again a measure of spurious currents and is shown in Figs. 4(a) and (b) for both formulations on the two different meshes. We observe that the spurious currents are larger for this case (as compared to imposing exact curvature) but remain bounded at large times when a balanced formulation is employed. This is in contrast to the conclusions in [17] where the use of CSF model with convolution of volume fraction for curvature calculation led to a constant increase in spurious currents even for a balanced formulation. The boundedness of the spurious currents in the present case, even with a less accurate curvature calculation is likely due to the specific approach adopted to effect the discrete balance. A comparison of the present algorithm with the results in [17] using convolved VoF is presented in Table 2 which shows the dependence of errors (L ∞ (u) and E ( ptotal )) on the timestep. It is clear that while the pressure errors and spurious currents are weakly dependent on timestep, the present results exhibit nearly one order lower spurious currents compared to the results of Francois et al. in [17]. The methodology in [17] which uses a collocated framework, involves interpolations of quantities from cells to faces (and vice-versa) while the present approach, on account of its

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

215

Table 2 The effect of timestep on spurious currents and pressure jump errors for inviscid static droplet with CSF convolution model on structured and unstructured meshes using balanced algorithm. The simulation is carried out till ρ t = 0.001 using ρh = 10 and W e = 1. l

t 10−3 10−4 10−5 10−6

Structured

Unstructured

Francois et al. [17]

L ∞ (u)

E ( ptotal )

L ∞ (u)

E ( ptotal )

L ∞ (u)

2.23 × 10−4 1.98 × 10−4 1.07 × 10−4 9.96 × 10−5

3.42 × 10−2 3.15 × 10−2 2.71 × 10−2 2.27 × 10−2

5.62 × 10−4 5.22 × 10−4 4.86 × 10−4 4.37 × 10−4

7.52 × 10−2 6.94 × 10−2 6.08 × 10−2 5.95 × 10−2

1.94 × 10−3 1.90 × 10−3 1.84 × 10−3 1.81 × 10−3

ρ

Fig. 5. Plots of the piezometric pressure superimposed with velocity vectors for ρh = 1000, in inviscid stationary water column with (a) unbalanced and l (b) balanced formulations. The scale for the vectors are identical in both figures.

hybrid staggered/non-staggered framework, achieves the balance in a relatively more straightforward manner as described in Section 4. Not surprisingly, the unbalanced formulation indeed shows large spurious currents that increase with time. These results highlight the importance of the balanced formulation and we believe that the present implementation would suffice to handle more complex problems, although improved curvature calculations using the approaches in [17,19] would definitely be beneficial in further lowering the magnitude of spurious currents. 6.2. Stationary water column in tank The importance of balanced formulation when body forces are dominant in multiphase flows is investigated using a stationary water column. A tank of unit height and width is filled up to half of its height with water. √ We assume an inviscid flow with a density ratio of 1000 and F r = 1. The length and velocity scales are chosen as W and g W , where W is the width of the tank. Simulations are carried out on a 100 × 100 Cartesian mesh using t = 0.001. The pressure contours at time t = 1 using the unbalanced and balanced formulations are shown in Figs. 5(a) and (b) respectively. It is evident that the unbalanced formulation generates significant spurious currents even at smaller times and predicts an erroneous pressure jump across the interface. The balanced formulation however estimates the exact pressure jump (which in this case is 499.5) and keeps the spurious currents to very low values of O(10−7 ) even at large times. The use of piezometric pressure, as also discussed earlier, leads to treatment of gravity as an interfacial force akin to surface tension. One can however achieve a balanced force approach even without adopting such a formulation. In the case a modified pressure is not employed, the body force would be a non-interfacial (volumetric) one. It would have a magnitude equal to ρ f g and would be lumped at the face centre which may be considered as the centroid of the momentum control volume as well (See Fig. 1). Following this approach, the pressure would be hydrostatic and not piezometric. This is reproduced in our numerical simulations as well, where Fig. 6 shows the correct linear pressure variation and no significant spurious currents are observed in this case even at t = 20. This is in sharp contrast to the conclusions in [21] where the cellcentred treatment of gravity, albeit in a collocated framework, resulted in significant spurious currents and pressure errors. This study therefore points to the fact that use of modified pressure and subsequently identical operators for pressure and body forces is not essential for a balanced force formulation. Nevertheless, we shall continue with the use of piezometric pressure and interfacial treatment of gravity in the studies to follow, unless otherwise specified. 6.3. Electromagnetic body force We further investigate the non-interfacial treatment of body force source terms by considering a discontinuous electromagnetic force. A 1 × 0.25 rectangular domain with a central region of length 0.3 assumed to be filled with ferro-fluid

216

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

Fig. 6. Pressure variation superimposed with velocity vectors for ρh = 1000, in inviscid stationary water column at t = 20. The velocities are O(10−7 ). l ρ

Fig. 7. (a) Computational domain with triangulated grid and (b) Pressure contours for a discontinuous electromagnetic body force. The velocities are of the O(10−8 ).

is discretised using a triangular mesh with 4010 elements (See Fig. 7(a)). A dimensionless electromagnetic force (per unit volume) of magnitude 100 is applied along the length in the region filled with ferro-fluid while the body force everywhere outside of this region is zero. The effects of gravity and surface tension are not included in this study. We define a linear distribution of density ρ = 1000 − 999(1 − x) in the domain and the simulation is carried out with a timestep of t = 0.001. The pressure variation, shown in Fig. 7(b) gives a total pressure jump of 29.99 while the analytical pressure jump is 30. In addition, very low spurious currents of the O(10−8 ) are observed, and these indicate that the stationary solution is accurately computed. It must be remarked that one could also define a modified pressure accounting for the electromagnetic body forces for this problem and then construct a well-balanced algorithm. This study reinforces the fact that a well-balanced algorithm does not necessarily require the use of modified pressures and employing identical discrete operators is not a pre-requisite for design of balanced force algorithms even on unstructured meshes. 6.4. Migration of a thermo-capillary droplet We demonstrate the efficacy of the well-balanced algorithm for multi-physics flow problems by considering the case of droplet migration under thermo-capillary effects. We consider a droplet of unit diameter centred at (2.5, 1.5) in a 5 × 7.5 rectangular domain which is subjected to a time-invariant linear temperature profile. The bottom and top walls are considered as no-slip boundaries and have dimensionless temperatures of 0 and 1 respectively while the lateral boundaries

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

217

Fig. 8. (a) Normalised rise velocity and (b) Shape of droplet for the thermo-capillary test case.

are slip surfaces. The surface tension is a function of temperature for this case and the source term due to surface tension forces now reads,

Fˇ st =

κ ∇φ We

−

1 Re

( T − T 0 )κ ∇φ −

1 Re

∇ T |∇φ| +

1 ∇ T · ∇φ Re

|∇φ|

∇φ

(22)

It must be noted that when implemented in the momentum equation, Eq. (13), the additional terms appearing in the surface tension force (which are due to thermal gradients) would still involve normal derivatives and these are then discretised using the same operator as described in Section 4. We choose Reynolds and Weber numbers of 0.06666 and 0.00444 respectively while density and viscosity ratios are taken as unity following [36]. A uniform grid with x = y = 0.02 is employed with a timestep of t = 0.001 and exact curvature is specified in this case. We do not consider effects of gravity in this study. Fig. 8(a) shows the rise velocity of the droplet normalised with the analytical value of 8.888 × 10−3 which is also the theoretical rise velocity of a spherical droplet for the zero Marangoni number limit in this study. Although we would expect the normalised rise velocity to be unity, the two-dimensional nature of the problem in this study gives a much smaller normalised terminal velocity of 0.825. This result is in close agreement with the RLSG scheme of [36] that reported a normalised rise velocity of 0.81 with a balanced force formulation on collocated grids. Furthermore, the history of rise velocity does not show any oscillatory behaviour for the VoF method as observed in [36] and settles to a constant terminal value in the present study. The droplet shape (superimposed with the velocity vectors) after terminal velocity has been attained is shown in Fig. 8(b) and the negligible deformation of the droplet which retains a circular shape justifies the imposition of exact curvature in this study. 6.5. Convection of large density droplet We now examine the importance of consistent formulation by considering the convection of a large density droplet in a stagnant fluid. A droplet of diameter D = 1 centred in a 3 × 3 domain is moved with constant unit velocity (U = 1) in an outer stagnant fluid. We choose a uniform grid with x = y = 0.01 and a density ratio of 106 . The timestep for the simulation is t = 10−4 and the simulation is carried out till t = 1. Gravity, surface tension and viscous effects are neglected in this study to focus solely on the discrete mass-momentum transport. The circular droplet shapes are significantly distorted in inconsistent formulation as shown in Fig. 9(a). This is due to dissimilar choice of convection schemes which leads to discretely inconsistent transport of mass and momentum and hence the spurious currents. The consistent formulation, on the other hand, largely preserves the shape of the droplet as the time progresses (shown in Fig. 9(b)) despite some diffusion of the interface which is likely due to the coarser mesh resolution. 6.6. Translating droplet Having discussed the effects of balanced formulation and consistent transport separately, we now discuss their combined effect in this subsection. A viscous drop is placed in a channel of length 5 and height 3 at (1, 1.5). The drop has a diameter D = 1 and the domain is discretised using a Cartesian mesh with x = y = 0.02. We specify a uniform inlet velocity of unity and Neumann boundary condition at outlet with the lateral surfaces being slip boundaries. The liquid inside the droplet is initialised to same velocity as the inlet and only surface tension effects are accounted for by specifying the exact curvature,2 whilst neglecting gravity.

2 It is not correct to apply the exact curvature when the drop deforms, but we still do so to avoid spurious errors that would arise from curvature calculations.

218

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

Fig. 9. Shape of the droplet of density ratio ρh = 106 moving at constant velocity U = 1 with timestep t = 10−4 using (a) inconsistent and (b) consistent l formulations. ρ

Studies are first performed at W e = 1 for two different Reynolds numbers and density ratios (10 and 1000) by considering fluids of equal dynamic viscosity. A timestep of t = 0.001 is chosen for the simulations which are carried out up to t = 3. Figs. 10(a)–(d) show the maximum error in the velocity in the domain relative to the background fluid flow and therefore is a direct measure of the spurious currents. We see from Figs. 10 (a) and (c) that for an inconsistent formulation, the spurious currents are nearly independent of whether a balanced approach is employed or not, for a given Reynolds number and density ratio. A more interesting fact evident from these figures is that although the spurious currents are the highest for the largest density ratio at Re = 10, they are quite significant even for a density ratio of 10 at Re = 1000. It is easy to see that the spurious currents in these two cases are equal or more than the background velocity itself, causing unphysical deformation of the droplet (not shown here for sake of brevity) as time progresses. On the other hand, the use of a consistent formulation exhibits much smaller spurious currents in comparison (as shown in Figs. 10(b) and (d)), with the use of a balanced force approach leading to a reduction in spurious currents by nearly an order of magnitude as compared to an unbalanced force approach. We also remark that only the use of a consistent formulation closely preserved the circular shape of the droplet in all cases. We conduct a second study to investigate the influence of viscosity ratio on the numerical results but keep the Reynolds number and density ratio constant at 10 and 1000 respectively. Figs. 11(a) and (b) show the results employing a balanced force approach with inconsistent and consistent transport respectively for viscosity ratios of 1 and 1000. We observe that the spurious currents are significantly smaller at larger (realistic) viscosity ratios, although the consistent formulation still produces lesser currents as compared to the inconsistent formulation. This may be possibly attributed to the pronounced damping at larger viscosity ratios. For the smaller viscosity ratio, as also observed in the earlier study, only the consistent approach preserves the shape of the droplet and results in acceptably low values of spurious currents. A third study is conducted at low Capillary numbers for static and moving droplets with both exact and calculated curvature. We choose W e = 0.01 and Re = 10, which corresponds to a capillary number of Ca = 0.001, with large density and viscosity ratios (both 1000). Simulations are performed using only the consistent and balanced formulation for the case of the static (specifying zero inlet velocity) and the translating droplet. The results in Figs. 12(a) and (b) show that the calculation of curvature introduces higher spurious currents as compared to exact specified interface curvature in both static and moving cases. Furthermore, the increase in spurious currents due to errors in curvature calculation are less pronounced in the moving case but are more prominent in the static case. While one could argue that the moving and static problems are essentially equivalent (differing by only the constant background velocity field), the calculated spurious errors appear to be significantly different, for a given curvature calculation. While this result is surprising and requires further investigation it appears to be consistent with the observations in the recent investigation of [15], who showed that the use of a moving reference frame in slug flows led to significantly reduced spurious currents as opposed to calculations on a stationary fixed

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

219

μ

Fig. 10. L ∞ norm of the velocity error relative to the background fluid moving with velocity U = 1 for W e = 1, μh = 1 and t = 0.001 using different l formulations: (a) inconsistent-balanced (b) consistent-balanced (c) inconsistent-unbalanced (d) consistent-unbalanced.

Fig. 11. L ∞ norm of the velocity error relative to the background fluid moving with velocity U = 1 for W e = 1, Re = 10 and t = 0.001 using (a) inconsistent-balanced and (b) consistent-balanced formulations.

reference frame with moving bubbles. Nevertheless, these studies clearly demonstrate that the consistent and balanced force algorithm performs the best leading to acceptably low spurious currents at high density ratios and even at low Capillary numbers where traditional algorithms are known to often experience difficulties. We also quantify the mass losses in the proposed algorithm by considering the translation of the droplet at W e = 1, Re = 10 with a density ratio equal to 1000 and unit dynamic viscosity ratio using the consistent transport and balanced force formulation. For this case we use a longer 12 × 3 computational domain with same grid resolution as employed in the previous studies. The upwind-biased CUI scheme is indeed prone to interface diffusion and we attempt to investigate this interface diffusion and associated mass losses in this numerical experiment. Figs. 13(a)–(f) show the interface diffusion

220

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

ρ

μ

Fig. 12. L ∞ velocity error of a droplet for W e = 0.01, Re = 10, ρh = 1000, μh = 1000 and t = 0.001 in a consistent-balanced formulation for (a) static l l and (b) moving cases.

ρ

Fig. 13. Line contours of volume fraction at time levels t = (a) 0 (b) 2 (c) 4 (d) 6 (e) 8 and (f) 10 moving with velocity U = 1 for Re = 10, ρh = 1000, l μ W e = 1, μh = 1 and t = 0.001 using the consistent-balanced formulation. l

at different time instants in the course of simulation. While one can clearly see that the interface diffuses over time, the diffusion is observed to be limited within 2–3 cells even at long times. We also calculate the mass errors as, N 

E mass =

i =1

|φcn c − φc0 c | N  i =1

(23)

φc0 c

where N is the total number of control volumes in the domain and superscripts n and 0 denote the numerical and initial volume fractions. The temporal history of E mass is presented in Fig. 14 and is acceptably low and remains bounded at long times. This study establishes the fact that the interface-free reconstruction approach in conjunction with consistent mass/momentum transport and well-balanced approach is an ideal methodology for computations of realistic flow phenomena with high density ratios.

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

ρ

221

μ

Fig. 14. Mass error of a translating droplet moving with velocity U = 1 for Re = 10, ρh = 1000, W e = 1, μh = 1 and t = 0.001 using the consistentl l balanced formulation.

Table 3 Simulation parameters for two fluid Poiseuille flow. Test

Re

ρh ρl

μh μl

Fr

We

t

C1 C2 C3 C4

10 10 10 1000

10 1000 1000 10

1 1 100 1

1 100 1 1

– – – –

0.001 0.001 0.001 1 × 10−6

6.7. Two fluid Poiseuille flow The Poiseuille flow of two immiscible fluids in a channel provides an ideal test case to highlight the importance of consistent formulation in gravity-dominated flows. We consider a horizontal channel with dimensions of [0, 5] × [−1, 1] which is discretised using a uniform Cartesian mesh with x = y = 0.02. We set the inlet and outlet pressures to zero √ and apply the gravity in the flow direction. The velocity scale is here chosen as g W and the width of the channel W is chosen as the length scale. We neglect the effects of surface tension in this study and carry out simulations for four different sets of parameters as shown in Table 3. The analytical fully developed profile may be easily derived as,

⎧  1 ⎨ Re + F r  2(1+β) U= 1 ⎩ Re + Fr 2(1+β) ρ

μ

2 βy −y 2(1+β)  2 y − 2(1+β)



α

−α 



1 2(1+β)

−



y if y ≥ 0 2(1+β)  y y2 − 2β if y < 0 2β(1+β)

1 2(1+β)

−

where α = ρh and β = μh represent the density and viscosity ratios respectively. l l The velocity profiles at the outlet obtained with both consistent and inconsistent formulations, which is balanced in either case, are shown in Figs. 15(a)–(d) for the four different cases. It is easy to see that the steady state solutions obtained using a consistent formulation agree well with the analytical result in all cases. The inconsistent formulation, on the other hand, performs as well as the consistent formulation in cases C1 and C3 but fails completely in cases C2 and C4. This is also evident in the significant distortion of the interface shown in Figs. 16(a) and (b). It must be remarked that while C2 represents a flow with high density ratio, case C4 involves a much smaller density ratio albeit at a larger Reynolds number. While cases C1, C2 and C4 involve unit viscosity ratio, Case C3 considers more realistic density and viscosity ratios and results in stable and accurate solution using both formulations, though at a lower Reynolds number. While the use of dissimilar schemes for convection lead to an inconsistent transport and unstable solutions for cases C2 and C4, one could argue that the different orders of accuracy of convective schemes used in the inconsistent formulation (1 for upwind and 2 for the CUI scheme) may be the source of instability. In order to investigate if this is indeed the case, we repeated the case C2 (which is more stringent) using dissimilar convection schemes of second-order accuracy for both mass and momentum transport. While the scheme N1 uses the high resolution scheme in [37] for momentum and CUI for volume fraction, scheme N2 adopts CUI for momentum and compressive CUIBS scheme for volume fraction [31]. The velocity profiles obtained using these schemes are presented in Fig. 17 and show that anomalous velocities are obtained in both cases, which underscore the fact that the dissimilar advection schemes themselves and not their orders of accuracy is responsible for the error. This also strongly suggests that the use of a consistent formulation may be necessary for high density multiphase flow simulations and in particular those involving high shear.

222

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

Fig. 15. Non-dimensional velocity profiles in two fluid Poiseuille flow for test cases (a) C1 (b) C2 (c) C3 and (d) C4.

Fig. 16. The distortion of the interface using consistent and inconsistent formulations for test cases (a) C2 and (b) C4.

6.8. Simulations of flow instabilities Flow instabilities are ubiquitous in many applications and typically involve high density ratios and/or high shear. We consider the simulations of Rayleigh–Taylor and Kelvin–Helmholtz instabilities to investigate the importance of consistent and inconsistent formulations in such flows. 6.8.1. Rayleigh–Taylor instability We choose a rectangular domain [−0.5, 0.5]√× [0, 4] for the simulations which is triangulated into 22256 elements. The velocity scale and length scale are chosen as g W and the width of the domain W respectively. In order to avoid any

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

223

Fig. 17. Non-dimensional velocity profiles in C2 for two fluid Poiseuille flow.

ρ

μ

Fig. 18. Time plots of Rayleigh–Taylor instability for ρH = 3, μH = 3, Re = 256 and F r = 1 on unstructured grid using (a) inconsistent and (b) consistent l l formulations.

effects of mesh asymmetry on the interface dynamics, we construct the unstructured mesh symmetric about the vertical centreline. The interface is initially perturbed as y = 2 + 0.1 cos(2π x) and time step is fixed to t = 0.001. We study the instability for two cases, of low and high Atwood numbers in moderate and strong gravity fields respectively. The first case following the numerical studies of Wang et al. [38] has an Atwood number 0.5 which corresponds to a density ratio of 3 and the viscosity ratio is also chosen as 3. The other dimensionless parameters for the simulation are Re = 256, F r = 1. The interface morphology at different time instants using inconsistent and consistent formulations are shown in the Figs. 18(a) and (b). The time evolution of the interfaces show excellent qualitative agreement while

224

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

Fig. 19. Time history of the interface position on left wall (top three lines) and at the vertical centreline (bottom three lines) for low At = 0.5.

ρ

μ

Fig. 20. Time plots of Rayleigh–Taylor instability for ρh = 7, μh = 1, Re = 1000 and F r = 0.1 on unstructured grid using (a) inconsistent and (b) consistent l l formulations.

remaining symmetric, independent of whether a consistent approach is employed. This is also reflected in the bubble and spike positions shown in Fig. 19, which agree well with the results of Wang et al. [38]. For the second case, we choose a strong gravity field corresponding to F r = 0.1 for a density ratio of 7, that translates to a moderately large Atwood number of 0.75. We choose a larger Re = 1000, while keeping the viscosity ratio equal to unity. The position of the interface at different instants of time using inconsistent and consistent formulations are shown in Figs. 20(a) and (b) respectively. It is evident that the inconsistent formulation converges to a non-physical solution with

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

ρ

225

μ

Fig. 21. The volume fraction profile of Kelvin–Helmholtz instability for ρh = 2, μh = 1 and Re = 5000 using (a) inconsistent and (b) consistent formulations. l l

a significant symmetry breaking at large times leading to a distorted interface. On the contrary, the consistent formulation preserves the symmetry about the vertical centreline correctly predicting a central spike even on the unstructured mesh. The results show that the use of an inconsistent formulation leads to a robust and accurate simulation when employed for a low Atwood number case, while it results in anomalous results for larger Atwood and Reynolds numbers in a stronger gravitational field. One can therefore conclude that the need for consistent formulations are not driven solely by the density ratio, but by a multitude of factors, including the viscosity ratio and the body forces. 6.8.2. Kelvin–Helmholtz instability We study the ability of the proposed numerical framework to deal with high shear flows by performing viscous simulations of shear instability in a unit square domain. The domain is discretised into a uniform 200 × 200 mesh and periodic boundary conditions are applied in the x-direction, while other boundaries are considered to be slip surfaces. The two fluids which have a density ratio of 2 and equal viscosities are separated by an interface at the mid plane corresponding to y = 0.5. We initialise the velocity field and volume fraction in the domain as follows [39].



u = tanh v =0

y − 0.5 − 0.01 sin(2π x)



φ = 0.5 1 + tanh



(24)

0.01



y − 0.5 − 0.01 sin(2π x) 0.01



(25) (26)

Simulations are carried out using t = 5 × 10−4 up to t = 1 using both consistent and inconsistent formulations. In these studies, we choose Re = 5000 and neglect the effects of gravity to focus on the effect of high shear. Figs. 21(a) and (b) show the interface between the fluids using the two approaches at t = 0.5 and t = 1 respectively. It is evident that the inconsistent formulation introduces errors which cause anomalous deformation of the interface while also leading to a completely different roll-up pattern, as compared to the consistent formulation.

226

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

Fig. 22. Comparison of the bubble shape at time t = 4.2.

ρ

μ

Fig. 23. Time evolution of (a) bubble mass centre (b) rising velocity for ρh = 1000, μh = 100, Re = 3.5, W e = 0.125 and F r = 1 on unstructured grid. l l

The results for the simulation of R–T and K–H instabilities further underline the importance of consistent massmomentum transport in flows with high density ratio as well as high shear. Coupled with a balanced force framework, we believe that a consistent formulation as described in this work, presents a robust approach for simulating complex multiphase flows involving large density contrasts and high shear such as those encountered in primary atomisation. 6.9. Rising bubble in stagnant fluid We now consider a multi-physics test case, involving both effects of gravity and surface tension, that has been the subject of many investigations using multiphase flow solvers. We simulate the rise of an air bubble in stagnant water using the benchmark studies of Hysing et al. [40] and choose Re = 3.5, W e = 0.125 and F r = 1 with the density and viscosity ratios being 1000 and 100 respectively. This test case is a challenging one for multiphase flow solvers, in particular because the bubble lies in a regime where break-up can possibly occur. The bubble of diameter D = 1 is initially positioned at (1, 1) in a computational domain of 2 × 4 which is discretised into 80426 triangular elements, ensuring grid symmetry √ about the centreline of the channel. The length and velocity scales are the diameter of the bubble D and g D respectively. Simulations are carried out up to t = 4.2 with a timestep of t = 0.001 with both consistent and inconsistent formulations. The rising bubble continuously deforms due to buoyancy and drag with the surface tension trying to provide a stabilising effect. The shape of the bubble at t = 4.2 shown in Fig. 22 is compared with the results (of TP2D code) in [40] and shows a very good agreement when employing either formulation. Interestingly, the evolution of mass centre of the bubble and rise velocity obtained from both formulations (See Figs. 23(a) and (b)) also agree well with the results in [40]. We now repeat this experiment with a viscosity ratio of unity while keeping all other parameters unchanged. The shapes of the bubble at three different time instants are shown in Figs. 24(a) and (b) for the inconsistent and consistent formulations respectively. While a validation study for this test is not available, the inconsistent formulation introduces non-physical distortion of the interface while the consistent formulation largely preserves the left-right symmetry of the bubble. The results from these numerical experiments show that the inconsistent formulations may lead to physically consistent solutions even for high density ratios under some conditions but could fail in other scenarios.

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

227

ρ

μ

Fig. 24. Shape of the rising bubble in (a) inconsistent and (b) consistent formulations at different times on unstructured grid for ρh = 1000, μh = 1, l l Re = 3.5, W e = 0.125 and F r = 1.

6.10. Droplet splashing on thin liquid film We now assess the importance of consistent transport in multiphase simulations by considering the problem of a droplet splashing on a thin layer of liquid. This problem is representative of real-life applications such as inkjet printing and involves complex flow physics, which necessitates a robust and accurate flow solver. For the simulations, we consider a 6 × 2 rectangular domain with a uniform grid resolution of x = y = 0.01. A thin layer of liquid is filled to one-tenth the domain height and a droplet of the same liquid with diameter D = 1, centrally placed just above the layer is given an initial downward velocity of U = 1 (downward). We consider that the working fluids are air and water which correspond to realistic density and viscosity ratios of 815 and 55 at room temperature. Since the impact occurs very quickly, gravity does not play a major role in this case and may be neglected. Simulations are carried out at Re = 66 and W e = 0.126 using a timestep t = 0.001 up to t = 1.5. The precess of splashing is pictorially represented in Figs. 25(a) and (b) respectively for inconsistent and consistent formulations. We observe significant differences between the two formulations with the former exhibiting unphysical asymmetry of the solution at later times (t > 1). The consistent formulation predicts a symmetric splash leading to the “crown” at t = 1.5 and the results agree favourably in the qualitative sense with those in Li et al. [41] (where the Re = 1000 based on heavier fluid properties). However, unlike in the results in [41] which were obtained using a Lattice Boltzmann method, we observe satellite droplets (which are diffused due to a relatively coarse mesh) and air entrapment in the thin layer. Fig. 26 shows the temporal evolution of the spread factor Dr , where r is the crown radius with the dimensionless

Ut . Despite being a two-dimensional simulation, the D U t 0.5 ( D ) , irrespective of the consistency of the formulation.

time

evolution of crown radius follows a power-law defined by

r D

=

We remark that this agreement for an inconsistent formulation is misleading, since it wrongly predicts the complex underlying physics of the impact. We further repeat this test case at a lower Re = 6.6 with all other conditions identical and results from both the formulations appear to agree reasonably well with one another and also with those reported in [41]. This can be observed in Figs. 27(a) and (b), which show the droplet splashing at the lower Re using the two formulations, where no secondary droplet formation is seen and a thicker liquid sheet emanates after impact. This corroborates the fact that while the consistent formulation can handle a wide variety of density ratios including high density contrasts, an inconsistent formulation could work equally well under a specified set of conditions. While not necessarily the case, these specified conditions can

228

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

ρ

μ

Fig. 25. Instantaneous profile of droplet splashing on thin liquid layer for ρh = 815, μh = 55, Re = 66 and W e = 0.126 using (a) inconsistent and (b) conl l sistent formulations.

correspond to some of the realistic scenarios (as in the bubble rising case, see Section 6.9) which might also explain why existing multiphase solvers have not deeply addressed the critical issue of consistent advection. 6.11. Collapse of a water column As a final test of balanced and consistent advection, we choose the widely researched problem of collapse of a water column [26,37], owing to its importance in free surface flows where large distortions of the interface are common. We discretise the domain which is a 7 × 2 rectangular tank using a triangulated mesh of 78668 elements. A 1 × 1 water column placed at the left corner surrounded by air, leading to density and viscosity ratios of 815 and 63 respectively. For the bottom wall we provide no-slip condition and the remaining tank walls are assumed to be slip boundaries. The Reynolds number for the simulation is chosen to be 2950, while Weber number is 0.54 and the Froude number is unity. The velocity √ and length scales in defining these dimensionless numbers are g H and the height of the water column H respectively. The simulation is carried out using t = 0.001 by employing consistent and inconsistent formulations up to a final time of t = 4.

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

229

Fig. 26. The evolution of spread factor with non-dimensional time for Re = 66.

ρ

μ

Fig. 27. Instantaneous profile of droplet splashing on thin liquid layer for ρh = 815, μh = 55, Re = 6.6 and W e = 0.126 using (a) inconsistent and (b) conl l sistent formulations.

230

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

ρ

μ

Fig. 28. Instantaneous profile of collapse of water column for ρh = 815, μh = 63, Re = 2950, W e = 0.54 and F r = 1 using (a) inconsistent and (b) consistent l l formulations.

Figs. 28(a) and 28(b) show the contour plots of density at time t = 0, 1, 2, 3, and 4 for inconsistent and consistent formulations respectively. While comparing both the formulations, it is seen that the density front (position of the interface at the bottom wall) using the inconsistent formulation moves at a slow rate in the initial stages and gets completely distorted in later stages. The consistent formulation, on the other hand, shows an error-free evolution of the front which is in good agreement with the experimental data of Martin and Moyce [42] as can be seen in Fig. 29. The time history of the evolution of the leading edge using the inconsistent formulation is obviously unphysical as can also be seen from Fig. 29 and is a direct consequence of dissimilar transport of mass and momentum at the discrete level. This study provides a strong justification for the need to employ consistent transport of mass and momentum in the simulation of high Re and high density ratio free surface flows which are of practical interest. 7. Consistent versus inconsistent transport for multiphase simulations The results from numerical investigations clearly show the need for a consistent transport of mass and momentum for accurate simulations of multiphase flows. Interestingly, a lack of consistent transport appears to yield accurate results for some cases while for other scenarios the numerical solutions are anomalous. It is therefore instructive to see under what conditions an inconsistent but balanced algorithm would suffice for stable and accurate simulations of multiphase flows. This is particularly relevant because many of the current multiphase flow solvers treat the convective terms in the LS/VoF equation and momentum equation with different schemes for discrete transport [19,23,40]. Table 4 summarises the results of various studies carried out in this work using the inconsistent formulation, with the relevant parameters, as discussed

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

231

Fig. 29. The front position of collapsing water column in air at the bottom of the tank.

Table 4 The simulation parameters for different test cases using the inconsistent formulation. Test case

ρh ρl

μh μl

Re

We

Fr

U

Physically consistent

Translating droplet

10 1000

1 1

10 10

1 1

– –

1 1

 ×

Poiseuille flow

10 1000 1000 10

1 1 100 1

10 10 10 1000

– – – –

1 100 1 1

– – – –

 ×  ×

Rayleigh–Taylor

3 7

3 1

256 1000

– –

1 0.1

– –

 ×

Kelvin–Helmholtz

2

1

5000

–

–

1

×

Bubble rising

1000 1000

100 1

3.5 3.5

0.125 0.125

1 1

– –

 ×

Droplet splashing

815 815

55 55

66 6.6

0.126 0.126

– –

1 1

× 

Collapse of water column

815

63

2950

0.54

1

–

×

in Section 6. In the last column  indicates a physically consistent result and one can see that the robustness of a solver incorporating an inconsistent approach is a function of the density and viscosity ratios as well as the non-dimensional numbers of interest which could include two or more of Re, W e and F r. It may be noticed that an inconsistent approach could work well for density ratios as high as 1000, as in one of the bubble rising test cases but could fail for comparatively low density ratios like 10, such as for the C4 case in two-fluid Poiseuille flow. A quick look reveals that the other parameters are however quite different in both these cases which suggests that difficulties in multiphase flow simulations do not stem merely from the presence of large density ratios. We can then conclude that numerical simulations with the inconsistent approach may or may not work satisfactorily at realistic density and viscosity contrasts, as witnessed for the bubble rise and droplet splashing test cases. The magnitude of spurious currents generated using an inconsistent transport is only known aposteriori and it is these numerically generated currents that decide the physical consistency of the numerical solution. On the contrary, the consistent approach along with a well-balanced algorithm guarantees the least spurious currents and fairly accurate non-anomalous solutions. In the light of extensive numerical experiments reported herein, the authors are of the opinion that a consistent and balanced force algorithm offers a robust and accurate approach for unsteady multiphase simulations, particularly those involving instabilities and complex physics. 8. Discussions and comparison with alternative algorithms The proposed algorithm is a novel approach for multiphase flow simulations that combines elements of balanced force treatment and consistent transport in a hybrid staggered/non-staggered framework. The algorithm is unique and different from other research efforts in the past and presents an advancement in terms of numerical techniques employed for multiphase flow simulations. We compare the proposed algorithm with other approaches in open literature for multiphase flows and highlight its salient advantages on three different counts.

232

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

1. On the balanced force implementation: Previous studies on balanced force approaches have been limited to surfacetension dominated flows on structured [17] and unstructured grids [19]. The importance of curvature calculations was emphasised in these studies to keep spurious currents to a minimum. While the efforts in [18] have looked at wellbalanced schemes in presence of body forces, these studies were restricted to structured meshes and also necessitating modifications in momentum interpolation. The proposed algorithm allows for a “natural” setting to enforce balancing and circumvents the use of momentum interpolation by virtue of the hybrid staggered/non-staggered framework. It is applicable to arbitrary polygonal meshes and can handle both interfacial and body forces without need for any special treatment pertaining to the forces. This is by virtue of the single “normal” momentum equation (Eq. (13)) where the pressure and other forces (such as those due to gravity/thermo-capillary effects) appear as normal derivatives. It must be noted that the Poisson equation for pressure correction (Eq. (15)) also features a normal derivative and the identical discretisation of this operator assures a discrete balance of forces with no further modifications. Moreover, we have also shown that a well-balanced algorithm for body forces can be realised even without the use of modified pressures and identical operators contrary to the work in [18]. The use of a hybrid staggered/non-staggered approach for multiphase flows therefore has the twin advantages of ease of applicability on unstructured meshes and simple construction of a well-balanced scheme. 2. On the consistent transport of mass and momentum: The need for using a consistent transport for mass and momentum has been recognised for long, but the specific implementation of the philosophy in numerical codes has been quite different in literature. While Rudman’s studies required the use of dual meshes [23], the approach of Raessi and Pitsch [26] was limited to staggered Cartesian meshes and employed a level-set methodology. The latter study defined a flux density which was used in momentum and mass transport equations and necessitated solution to an additional transport equation for density at cell faces. The extension of this approach to three-dimensional unstructured meshes is clearly a compute-intensive and non-trivial effort, which could limits its applicability to realistic scenarios. The proposed algorithm, which is motivated by the work of Ghods and Herrmann [20], overcomes these implementational hurdles while providing for a robust and accurate approach for simulating high density ratio flows. Similar to [20], we employ the same identical scheme for evaluating the convective fluxes in the mass and momentum equations, thereby incurring small but consistent errors in mass and momentum conservation. We however use a volume-of-fluid approach devoid of interface reconstruction along with bounded upwind-biased convective schemes while the studies in [20] employed first-order upwind schemes in conjunction with a level set methodology. The use of identical, nominally second-order accurate and monotone schemes for advection with an algebraic VoF method ensure accurate and robust simulations of multiphase flows with high density ratio and/or high shear. It must also be remarked that no special treatments such as the density-weighting in [19] are necessary for simulating large density ratio flows using the proposed algorithm. 3. On utility and generic approach: The studies of Francois et al. [17], Herrmann [43] as well as Denner and van Wachem [19] are focused on developing balanced-force approaches for surface-tension dominated multiphase flows, with the latter two studies applicable to unstructured meshes as well. While the efforts of Montazeri and Ward [18] tackle other forces as well, these are restricted to Cartesian meshes only. Stable simulations of high density ratio flows has been the subject of investigation in the studies of Raessi and Pitsch [26] as well as Ghods and Herrmann [20], with the former restricted to structured meshes. The present study considers well-balanced schemes as well as robust simulations of large density flows, both in isolation and in combination with an emphasis on realistic scenarios. Moreover, the proposed algorithm is equally applicable to any grid topology and is generic enough to accommodate interfacial and non-interfacial forces. The extensive numerical studies presented in this work which encompass a large number of problems including canonical cases and realistic scenarios demonstrate the generic nature of the algorithm and its viability for practical applications involving complex flows. The key features that make the present algorithm unique for multiphase flow simulations are the specific implementations of force balancing and consistent transport. The novelty of the hybrid staggered/non-staggered framework also facilitates a simple implementation of balanced-force approaches which can also simulate flows with large density ratios and/or shear when identical high-resolution schemes are used for convective transport. The studies show that the proposed approach is an elegant and “fail-safe” implementation for flows with high shear and/or density contrasts for realistic simulations of complex multiphase flows which maybe dominated by gravity, surface tension and/or thermo-capillary effects. 9. Summary and conclusions A novel and generic algorithm to simulate multiphase flows with large density ratios in the presence of dominant interfacial and/or body forces is proposed. The proposed algorithm employs a novel staggered/non-staggered framework, amenable on arbitrary polygonal meshes, which allows for a natural implementation of the balanced force framework. An accurate and consistent transport of mass and momentum is also realised by employing identical high resolution schemes for convective transport of these quantities at the discrete level. The algorithm which also employs an interface reconstruction-free VoF approach is tested rigorously on several multiphase flows of varying complexity that involve both gravity and surface tension dominated fluid flows using structured and unstructured grids with different density and viscosity contrasts. The salient conclusions from the numerical investigations may be summarised as follows.

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

233

1. The use of identical discrete operators is not always necessary to ensure a well-balanced algorithm, particularly when body forces such as gravity are involved. Our studies have shown that gravity may be as an interfacial force (if piezometric pressure is employed) or as a non-interfacial volumetric force (if hydrostatic pressure is employed). While the former necessarily requires identical discrete operators for gradients of pressure and density, the latter does not require any such treatment. The present study conclusively demonstrates that a non-piezometric approach to handle gravitational forces (and other similar non-interfacial forces) can also result in a balanced-force approach which is in contrast to existing studies. 2. A consistent discretisation of convective terms in volume fraction and momentum equations is important to ensure that the simulations for high density ratio flows remain stable at all times and the numerical solutions are physically consistent. This can be realised by using identical convection schemes with the high-resolution schemes in the present study providing nominal second-order accuracy in addition to robustness. 3. The use of balanced-force approach and consistent mass and momentum transport is desirable for realistic simulations of multiphase flows. However, for most cases the parasitic (spurious) currents from the lack of consistent transport is larger in magnitude when compared with those due to unbalanced approaches. While parasitic currents could introduce errors in the numerical solution, those due to a lack of consistent discrete transport can be detrimental to numerical stability as well. 4. The success of an inconsistent approach to mass and momentum transport is a complex function of density and viscosity ratios in addition to the Reynolds, Froude and Weber numbers. We speculate that this happens because the numerical diffusion suffices to keep the spurious currents bounded at a value lower than the bulk fluid velocity. While such an approach could work for some scenarios, it fails miserably in others and is therefore not generally recommended for multiphase applications. 5. The consistent transport coupled with a balanced force algorithm provides a numerically stable and reasonably accurate methodology to simulate realistic multiphase flows. The lack of interface reconstruction is an added advantage which makes extension of the proposed algorithm to three-dimensional unstructured grids relatively straightforward at a low computational cost. The authors believe that the proposed algorithm and computational framework offers a promising platform for simulations of realistic multiphase flows of interest. Moreover, while the present study considers an unique hybrid staggered/nonstaggered framework, it is possible to adopt the ideas propounded herein to collocated grids as well. The development of a fully collocated multiphase flow solver rooted in a scientific philosophy similar to that described in the present work is currently in progress and its application for multiphase flows shall be reported in a separate publication. Extension of the proposed framework to three-dimensions and simulation of more complex flows will also be undertaken in the near future. Appendix A. Derivation of the normal momentum equation and fractional step approach We provide here a concise derivation of the normal momentum equation and description of the fractional step approach in the hybrid staggered/non-staggered framework. Using Eq. (10) in Eq. (4) and representing the convective and viscous flux terms by the vectors Fc and F v respectively gives,

∇ρ ∂(ρ u) κ ∇φ ˇ + ∇ · Fc = −∇ p + ∇ · F v − + Fb e·x+ ∂t Fr We

(A.1)

The above vector conservation law is first projected along the direction of the normal to the face and then integrated over the control volume formed by the union of two cells sharing that face (See Fig. 1). The temporal term then becomes,



d(ρ U f ) ∂(ρ u) · n f d ≈  ∂t dt

(A.2)



where U f = u · n f is the normal velocity at the face and we assume that the time derivative remains invariant over the control volume. The temporal derivative is discretised using a three-point backward differencing scheme. The convection term which is in divergence form becomes,

 

   Fc · n d · n f ∇ · Fc · n f d ≈



(A.3)



where the Gauss divergence theorem has been applied and  denotes the control surface of the momentum control volume. Using the definition of the convective flux and using a single-point Gauss quadrature to evaluate the integral gives,

 



Fc · n d · n f =





ρe ue U e se · n f

(A.4)

e

where e represents the edges of the momentum control volume, which is union of two cells (See Fig. 1), while f denotes the face shared by the two cells. The term in parentheses is first evaluated with ue obtained from the upwind-biased

234

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

scheme as would be done in a collocated framework and then the dot product of the evaluated vector with n f is carried out resulting in the scalar contribution denoted as CONV in Eq. (13). An analogous approach for the viscous terms gives the diffusive contribution denoted as VISC in Eq. (13). We then have for the viscous term,

 

   F v · n d · n f ∇ · F v · n f d ≈



(A.5)



where F v =



2μ Re

D where D = 12 (∇ u + ∇ uT ) is the deformation tensor. Therefore we have,



F v · n d · n f =

  2μ

e

e



Re

 De · ne se · n f

(A.6)

where the deformation tensor is evaluated at the faces of the momentum control volume using a central differencing approach. It must be remarked that the deformation tensor would involve gradients of velocities and these are obtained using the Green–Gauss reconstruction [34]. Notice that the viscous flux computation is also done in the same manner as in non-staggered meshes. The pressure gradient term may be simplified as,

 ∇ p · n f d ≈ ∇ p f · n f  ≈

δ p    δn f

(A.7)



where we assume that the pressure gradient at the face is an averaged value over the momentum control volume and the normal pressure derivative is evaluated at the face centre. The surface tension term becomes,



κ ∇φ We

· n f d ≈

κf We

∇φ f · n f  ≈

κ f δφ  W e δn

 

(A.8)

f



We assume here that ∇φ f represents an averaged value over the control volume which is lumped at the face centre. Moreover, any similar additional terms (such as in case of thermo-capillary effects) are treated in a similar manner. The gravitational effects are restricted to the interface because of the use of a piezometric pressure and can be simplified as,



∇ρ Fr

(e · x) · n f d ≈

1 Fr

(e · x) f ∇ ρ f · n f  ≈

1 Fr

(e · x f )

δ ρ    δn f

(A.9)



where we again assume that the face-centred gradient is representative of a control-volume averaged quantity. The body force terms (other than gravity) are easily simplified as,



Fˇ b · n f d ≈ Fˇ bf · n f 

(A.10)



where the body force evaluated at the face centre is an averaged value over the momentum control volume. Note that the centre of face f is assumed to be the centroid of the momentum control volume in this derivation. All the above simplifications put together gives the scalar normal momentum equation on unstructured meshes which is described in a fully discrete form by Eq. (13), where the ‘ ’ term denotes the auxiliary normal momentum. The fractional step approach in this work is similar to that employed for single-phase incompressible fluid flows, except that we consider spatially varying density. The difference between the normal momentum equations at (n + 1) and intermediate (‘ ’) time levels is two-fold. Firstly, the normal velocity and centroidal velocity components are at the intermediate levels in the auxiliary momentum equation. Secondly, the pressure derivative is at time level n in the predictor step which gives the auxiliary momentum. Subtracting the discrete normal momentum equations and neglecting differences due to convective and viscous contributions between (n + 1) and ‘ ’ time levels gives,

3 (ρ U f )n+1 − 3 (ρ U f )

2 t

=−

δ( pn+1 − pn ) δn

ρ = ρ n+1 and dividing by the density gives,   U nf +1 − U f 2 1 δP =− t 3 ρ n +1 δ n

(A.11)

Realising that

(A.12)

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

235

Fig. 30. The shape of the Zalesak’s disk at the initial time (dashed line) and after one resolution (solid line) on unstructured mesh.

Fig. 31. Shape errors of a Zalesak’s disk after one revolution for different time steps. The dashed line represents the second order accuracy.

A discrete summation over the faces of a cell gives,



U f s f =

f

2 t  3

f

1

δP

ρ n +1 δ n

s f

(A.13)

which is the discrete Poisson equation for pressure correction given by Eq. (15). Notice that we have used the relation  f U f s f = 0 which is the discrete analogue of the incompressibility constraint in Eq. (1). We also emphasise that the control volumes for scalar quantities such as the pressure correction solved for by the Poisson equation are the individual cells themselves. The solution to the non-linear auxiliary momentum equation and linear pressure correction equation give the intermediate normal momentum and pressure corrections respectively. The pressure and normal momentum at the time level (n + 1) follow from algebraic corrections which are given by Eqs. (16) and (17). Appendix B. Interface capture and Zalesak disk problem We study the efficacy of the interface capturing scheme as applied to the volume fraction advection equation by considering the Zalesak solid body rotation problem [44]. This test is carried out by solving the advection equation, Eq. (18), with a prescribed velocity field so that it is decoupled from the Navier–Stokes equations. A disk of unit diameter is cut by a slot of 0.12 unit width and is centred at (2,2.25) in a 4 × 4 domain which is discretised into 40532 triangular elements. The disk is rotated with the velocity (u , v ) = (−0.5( y − y 0 ), 0.5(x − x0 )) about (2, 2) and the shape of the disk at t = 0 and after one revolution are plotted in Fig. 30. While the interface is subjected to some diffusion, expected of an upwind-biased scheme for advection, the shape remains largely symmetric. We also study the shape errors on this grid for four different time steps viz. t = 0.005, 0.0025, 0.00125 and 0.000625 after one revolution. The shape error is defined as, E shape =

N 

i =1

|φcn c − φce c |/

N 

i =1

φc0 c where superscripts n and e represent the numerical and exact solutions at the end of

one revolution while ‘0’ denotes the solution at initial time. It can be seen from Fig. 31 that the shape errors decrease at a rate of 1.59 as time step becomes smaller, indicating a temporal accuracy between 1 and 2 (but closer to second order) for the algebraic volume-of-fluid method in this work.

236

J.K. Patel, G. Natarajan / Journal of Computational Physics 350 (2017) 207–236

References [1] A. Prosperetti, G. Tryggvason, Computational Methods for Multiphase Flow, Cambridge University Press, 2009. [2] C.W. Hirt, B.D. Nichols, Volume of fluid (vof) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1) (1981) 201–225. [3] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1) (1988) 12–49. [4] S. Van der Pijl, A. Segal, C. Vuik, P. Wesseling, A mass-conserving level-set method for modelling of multi-phase flows, Int. J. Numer. Methods Fluids 47 (4) (2005) 339–361. [5] M. Sussman, E.G. Puckett, A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows, J. Comput. Phys. 162 (2) (2000) 301–337. [6] T. Wang, H. Li, Y. Feng, D. Shi, A coupled volume-of-fluid and level set (voset) method on dynamically adaptive quadtree grids, Int. J. Heat Mass Transf. 67 (2013) 70–73. [7] P.R. Redapangu, S. Vanka, K.C. Sahu, Multiphase lattice Boltzmann simulations of buoyancy-induced flow of two immiscible fluids with different viscosities, Eur. J. Mech. B, Fluids 34 (2012) 105–114. [8] J. Kim, A diffuse-interface model for axisymmetric immiscible two-phase flow, Appl. Math. Comput. 160 (2) (2005) 589–606. [9] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, G. Zanetti, Modelling merging and fragmentation in multiphase flows with surfer, J. Comput. Phys. 113 (1) (1994) 134–147. [10] D. Gerlach, G. Tomar, G. Biswas, F. Durst, Comparison of volume-of-fluid methods for surface tension-dominant two-phase flows, Int. J. Heat Mass Transf. 49 (3) (2006) 740–754. [11] S.J. Cummins, M.M. Francois, D.B. Kothe, Estimating curvature from volume fractions, Comput. Struct. 83 (6) (2005) 425–434. [12] S. Popinet, S. Zaleski, A front-tracking algorithm for accurate representation of surface tension, Int. J. Numer. Methods Fluids 30 (6) (1999) 775–793. [13] T. Abadie, J. Aubin, D. Legendre, On the combined effects of surface tension force calculation and interface advection on spurious currents within volume of fluid and level set frameworks, J. Comput. Phys. 297 (2015) 611–636. [14] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, J. Comput. Phys. 228 (16) (2009) 5838–5866. [15] Z. Pan, J.A. Weibel, S.V. Garimella, Spurious current suppression in vof-csf simulation of slug flow through small channels, Numer. Heat Transf., Part A, Appl. 67 (1) (2015) 1–12. [16] Y. Renardy, M. Renardy, PROST: a parabolic reconstruction of surface tension for the volume-of-fluid method, J. Comput. Phys. 183 (2) (2002) 400–421. [17] M.M. Francois, S.J. Cummins, E.D. Dendy, D.B. Kothe, J.M. Sicilian, M.W. Williams, A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework, J. Comput. Phys. 213 (1) (2006) 141–173. [18] H. Montazeri, C. Ward, A balanced-force algorithm for two-phase flows, J. Comput. Phys. 257 (2014) 645–669. [19] F. Denner, B.G. van Wachem, Fully-coupled balanced-force vof framework for arbitrary meshes with least-squares curvature evaluation from volume fractions, Numer. Heat Transf., Part B, Fundam. 65 (3) (2014) 218–255. [20] S. Ghods, M. Herrmann, A consistent rescaled momentum transport method for simulating large density ratio incompressible multiphase flows using level set methods, Phys. Scr. 2013 (T155) (2013) 014050. [21] H. Montazeri, M. Bussmann, J. Mostaghimi, Accurate implementation of forcing terms for two-phase flows into simple algorithm, Int. J. Multiph. Flow 45 (2012) 40–52. [22] J. Mencinger, I. Žun, On the finite volume discretization of discontinuous body force field on collocated grid: application to vof method, J. Comput. Phys. 221 (2) (2007) 524–538. [23] M. Rudman, A volume-tracking method for incompressible multifluid flows with large density variations, Int. J. Numer. Methods Fluids 28 (2) (1998) 357–378. [24] O. Desjardins, V. Moureau, Methods for multiphase flows with high density ratio, in: Center for Turbulent Research, Summer Programm 2010, 2010, pp. 313–322. [25] M. Bussmann, D.B. Kothe, J.M. Sicilian, Modeling high density ratio incompressible interfacial flows, in: ASME 2002 Joint US–European Fluids Engineering Division Conference, American Society of Mechanical Engineers, 2002, pp. 707–713. [26] M. Raessi, H. Pitsch, Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method, Comput. Fluids 63 (2012) 70–81. [27] J. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (2) (1992) 335–354. [28] O. Ubbink, Numerical Prediction of Two Fluid Systems with Sharp Interfaces, Ph.D. thesis, University of London, 1997. [29] L. Ge, F. Sotiropoulos, A numerical method for solving the 3d unsteady incompressible Navier–Stokes equations in curvilinear domains with complex immersed boundaries, J. Comput. Phys. 225 (2) (2007) 1782–1809. [30] G. Natarajan, F. Sotiropoulos, Adaptive finite volume incompressible Navier–Stokes solver for 3d flows with complex immersed boundaries, in: APS Division of Fluid Dynamics Meeting Abstracts, vol. 1, 2009. [31] J.K. Patel, G. Natarajan, A generic framework for design of interface capturing schemes for multi-fluid flows, Comput. Fluids 106 (2015) 108–118. [32] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, V. Eijkhout, W. Gropp, D. Kaushik, M. Knepley, et al., Petsc users manual revision 3.5, Argonne National Laboratory (ANL). [33] Lis libraries, http://www.ssisc.org/lis, 2013. [34] J. Blazek, Computational Fluid Dynamics: Principles and Applications, vol. 1, Elsevier, 2001. [35] N.P. Waterson, H. Deconinck, Design principles for bounded higher-order convection schemes-a unified approach, J. Comput. Phys. 224 (1) (2007) 182–207. [36] M. Herrmann, J. Lopez, P. Brady, M. Raessi, Thermocapillary motion of deformable drops and bubbles, in: Proceedings of the Summer Program, 2008, p. 155. [37] Y.-Y. Tsui, S.-W. Lin, T.-T. Cheng, T.-C. Wu, Flux-blending schemes for interface capture in two-fluid flows, Int. J. Heat Mass Transf. 52 (23) (2009) 5547–5556. [38] Y. Wang, C. Shu, H. Huang, C. Teo, Multiphase lattice Boltzmann flux solver for incompressible multiphase flows with large density ratio, J. Comput. Phys. 280 (2015) 404–423. [39] H.G. Lee, J. Kim, Two-dimensional Kelvin–Helmholtz instabilities of multi-component fluids, Eur. J. Mech. B, Fluids 49 (2015) 77–88. [40] S.-R. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, L. Tobiska, Quantitative benchmark computations of two-dimensional bubble dynamics, Int. J. Numer. Methods Fluids 60 (11) (2009) 1259–1288. [41] Q. Li, K. Luo, X. Li, Lattice Boltzmann modeling of multiphase flows at large density ratio with an improved pseudopotential model, Phys. Rev. E 87 (5) (2013) 053301. [42] J. Martin, W. Moyce, Part IV. An experimental study of the collapse of liquid columns on a rigid horizontal plane, Philos. Trans. R. Soc. Lond. Ser. A, Math. Phys. Eng. Sci. 244 (882) (1952) 312–324. [43] M. Herrmann, A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids, J. Comput. Phys. 227 (4) (2008) 2674–2706. [44] M. Rudman, Volume-tracking methods for interfacial flow calculations, Int. J. Numer. Methods Fluids 24 (7) (1997) 671–691.