A generic framework for design of interface capturing schemes for multi-fluid flows

Computers & Fluids 106 (2015) 108–118

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

A generic framework for design of interface capturing schemes for multi-fluid flows Jitendra Kumar Patel, Ganesh Natarajan ⇑ Department of Mechanical Engineering, Indian Institute of Technology Guwahati, India

a r t i c l e

i n f o

Article history: Received 20 April 2014 Received in revised form 24 August 2014 Accepted 3 October 2014 Available online 22 October 2014 Keywords: Interface capturing Volume-Of-Fluid Unstructured grids High resolution Flux limiters Boundedness

a b s t r a c t A generic methodology for the design of interface capturing schemes for immiscible multi-fluid flows using the VOF (Volume-Of-Fluid) approach is outlined. Interface capturing schemes are devised as a blend of high-resolution and compressive schemes and their efficiency and accuracy are dependent on the choice of the constituent schemes and the blending function. On the basis of a set of design principles proposed in this work, we show that interface capturing schemes proposed in literature may be encompassed into a single class of GPL (Generalised Piecewise Linear)-j schemes allowing for a unified approach for development of such schemes. This methodology is used to de-construct four existing schemes and to propose two new schemes for interface capture in an unstructured finite volume framework. The schemes are tested on challenging advection problems using both structured and unstructured grids to evaluate their performance. The new schemes perform as well as existing schemes and even outperform them on unstructured grids, while also exhibiting near Courant number independence for all test problems and grid topologies. The CUIBS (Cubic Upwind Interpolation based Blending Scheme) proposed in the study shows the best performance among all interface capturing schemes discussed herein and is applied to the study of immiscible binary fluid flows. Numerical simulations of viscous sloshing in a tank and Rayleigh–Taylor instability demonstrate the ability of the CUIBS scheme to resolve fluid–fluid interfaces with minimal diffusion. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction The success of numerical simulations of incompressible immiscible fluid flows depends on the ability of the algorithm to resolve the fluid interfaces accurately and preserve their sharpness. The Volume-Of-Fluid (VOF) method [1] is among the approaches employed for such simulations owing to its robustness and conservation properties. This approach involves solving advection equation(s) for the volume fraction(s) that define the interface(s) and utmost care needs to be taken to design stable and robust convective schemes that result in minimal numerical diffusion. While effective geometric methods to locally reconstruct the interface have been considered in literature [2,3], their implementation on unstructured meshes is a major challenge. The need for advection schemes that preserve a sharp interface without compromising on solution accuracy and allows for easy implementation on arbitrary mesh topologies leads to the concepts of interface capturing ⇑ Corresponding author at: D-308, Department of Mechanical Engineering, Indian Institute of Technology Guwahati, Guwahati 781039, Assam, India. Tel.: +91 361 2582685. E-mail address: [email protected] (G. Natarajan). http://dx.doi.org/10.1016/j.compfluid.2014.10.005 0045-7930/Ó 2014 Elsevier Ltd. All rights reserved.

schemes. Interface capturing schemes are constructed by combining high-resolution schemes with compressive schemes, using a suitable blending function, thereby obviating the need for any geometric reconstruction of the interfaces. Among the earliest works in this regard are the HRIC scheme of Muzaferija et al. [4] and the CICSAM scheme proposed by Ubbink and Issa [5]. A prominent drawback of both these methods is the deterioration of their performance with Courant number necessitating the use of smaller values of Courant number. Darwish and Moukalled [6] developed the STACS scheme to successfully overcome this drawback and subsequently proposed a family of transient interface capturing schemes referred to as TICS family of schemes [7] that combine a bounded high-order transient scheme with a bounded high-order compressive scheme. Gopala and co-workers [8] proposed a scheme akin to the Inter-Gamma scheme of Jasak and Weller [9] and investigated its utility for free-surface flows. Walters and Wolgemuth [10] proposed the Bounded Gradient Maximization scheme which involves a cell-based flux limiting with a suitable weighting factor while Xiao and co-workers [11] developed the THINC scheme based on the tangent function to compute volume fraction flux that ensures oscillation-free solution for fluid interface problems. Notable recent attempts in this direction are the

109

J.K. Patel, G. Natarajan / Computers & Fluids 106 (2015) 108–118

FBICS scheme [12] and the HiRAC scheme [13], with the latter using a modified advection equation with an artificial compression term to ensure sharp resolution with stability. Most of these schemes have been applied to time-dependent interfacial flows and have proven successful in resolving the complex flow dynamics. It must be remarked that there exist other approaches for multi-fluid and multiphase simulations in literature, most notably the Level Set (LS) approach by Osher and Sethian [14] where the interface is defined by a zero distance function and immiscible fluids on either side of interface are distinguished by positive and negative values of the distance function. While many researchers have successfully employed the LS approach in solving a wide variety of problems [15], its applications have mostly been restricted to structured Cartesian meshes and extension to unstructured meshes is not straightforward. An interesting approach to solving immiscible flows that has received significant attention is the diffuse interface technique [18], where the Navier– Stokes equations are solved in conjunction with convective Cahn– Hilliard equation. The method is conservative in nature but the interface needs special treatment to maintain its sharpness. Lattice Boltzmann Methods (LBM) have also been used for efficient simulation of multi-fluid flows which include studies of Rayleigh–Taylor instability in titled channel [16] and viscosity stratified flows [17]. The focus of the present work is on the VOF method with interface capturing as a methodology to handle multi-fluid immiscible flows on unstructured meshes in two dimensions. We emphasise that while several interface capturing schemes have been proposed in the past and are in use, there is no definitive framework to analyse these schemes for possible improvements and development of newer schemes. While almost all interface capture schemes are constructed as a blend of high-resolution and compressive schemes, there are no clear guidelines for the specific choices in open literature. Inspired by the unified approach for convective schemes introduced in [19] which brought several different higher-order convection schemes under a common framework, the authors have attempted to devise a comprehensive set of design principles for evaluation and development of interface capturing schemes for multi-fluid flows in this work. The rest of the paper is organised as follows. Section 2 briefly discusses the VOF approach and solution of the interface advection equation. Highresolution schemes and boundedness are discussed in Section 3 and the design principles for interface capturing schemes are proposed in Section 4. Examination of existing interface capture schemes and development of novel alternatives constitute Sections 5 and 6 and their performance for typical advection problems is the focus of Section 7. Section 8 describes the use of new interface capturing schemes for practical engineering problems such as liquid sloshing and Rayleigh–Taylor instability. A summary of the contributions of the study and possible future directions are given in Section 9. 2. Interface advection equation The Volume of Fluid (VOF) approach captures the interface by solving a transport equation for the scalar volume fraction. The interface between any two fluids which are immiscible will convect as a passive scalar and therefore satisfies a hyperbolic equation that reads,

solves (k  1) advection equations for as many volume fractions where the velocity field is determined from the solution to momentum equations. In this study, we restrict ourselves to binary fluid flows (k = 2) which are immiscible as well as incompressible. Consequently, we now solve for a single volume fraction and the discrete counterpart of the advection equation in an unstructured finite volume framework may be written in conservative form as,

Xc

X nþ1 n 3/nþ1  4/nc þ /n1 c c þ /f U f Dsf ¼ 0 2 Dt f 2C

where Xc is the volume (or area in two-dimensions) of the cell C; /c is the cell-averaged value of the volume fraction in that cell, /f denotes the value at the midpoint of faces constituting the cell, U f is the normal velocity at the faces constituting the cell and Dsf is face area. The value of the volume fraction obtained in each cell is limited to lie in the range [0, 1]. It must be remarked that P2 k k¼1 /c ¼ 1 and that the divergence-free condition is employed to recast the advection equation in its conservative form. The choice of second-order accurate backward differencing scheme (where the superscript now denotes the time level) is chosen to achieve time-accurate solutions for unsteady flows that are routinely encountered with immiscible binary fluids. The advection equation is solved assuming a constant value for the velocity field at a particular timestep, and is known from the solution to the momentum equations in the previous time-step. This approach of freezing the velocity field linearises the interface advection equation and the resulting system of linear algebraic equations (at a discrete level) are solved implicitly using an ILU-preconditioned GMRES approach [20]. It must be however remarked that for the tests carried out in Section 7, the (analytical) velocity field is prescribed directly, rather than computed as the solution to the momentum equations. 3. High resolution schemes and convection boundedness criterion The solution to the interface advection equation necessitates a robust convection scheme that ensures sharp resolution of the discontinuity. First order schemes are robust but they are overly diffusive and smear the interfaces. While higher-order (secondorder or more) accurate schemes appear to be a natural solution, they suffer from loss in monotonicity as a consequence of the Godunov Barrier Theorem [21]. The development of high-resolution schemes, which are at least second-order accurate with solution boundedness, has been a subject of research over the years. These schemes overcome the barrier theorem by introducing non-linear ‘‘flux limiters’’ to achieve solution monotonicity. Thus, the face value /f of the volume fraction on the uni-dimensional grid shown in Fig. 1(a) may be determined as,

/f ¼ /C þ

/C  /U  cðrÞ 2

where /C and /U are values of volume fraction at cell centroid of the upwind and far upwind cell respectively. Here, cðrÞ represents the flux limiter (FL) which depends on the gradient ratio r, which

U rCD

k

@/ þ u  r/k ¼ 0 @t k

ð1Þ

The quantity / represents the fraction of the kth fluid contained in the cell, where the fluids are arranged in increasing order of densities. Therefore, for the case of k immiscible fluids, this approach

ð2Þ

C U

c

C

(a)

f

D

f

D

U (b)

Fig. 1. Grid topology (a) structured mesh and (b) unstructured mesh.

110

J.K. Patel, G. Natarajan / Computers & Fluids 106 (2015) 108–118

may be defined as the ratio of centred to upwind gradient. In order to have non-oscillatory and robust convection schemes, the flux limiter must be chosen so as to ensure boundedness. The Convection Boundedness Criterion (CBC) [22] which is based on normalised variable approach (NV) is chosen in this study. This criterion is quite robust in practice and amenable for application on arbitrary polygonal meshes such as those encountered in this study. The normalised value of the volume fraction is defined by,

~ ¼ /  /U / /D  /U Recall that the gradient ratio r is defined by,

r¼

~f ¼ / ~ f ðBDÞ f ðhf Þ þ / ~ f ðHRÞ ð1  f ðhf ÞÞ /

/D  /C /C  /U

where

It is therefore easy to show that the FL and NV approaches are related by,

~C ¼ /

1 rþ1

1. Blending function: The blending function must be a continuous function that allows for a smooth transition between the HR and BD schemes which is dictated by the orientation of the physical interface with the mesh. A sharp interface will result when the physical interface is aligned with the cell face if the fully compressive scheme is used, while a purely high-resolution scheme is desirable when the cell face is normal to the physical interface. Of the different blending functions proposed in previous works, the function f ðhf Þ ¼ cos4 hf is found to work well and is adopted in the present study. A computationally less expensive blending devoid of trigonometric functions as suggested in [13] may also be pursued. Therefore, the normalised face value for the volume fraction may be obtained as,

  ~f ¼ / ~ C 1 þ cðrÞ / 2

ð3Þ

The FL/NV transformations described above is particularly helpful to understand the construction of interface capturing schemes as discussed in the later sections. The CBC criterion is based on the idea that the normalised face values be bounded by the normalised values in cells sharing the face, which in turn implies that the ~ f ¼ 1, boundedness region is a triangle enclosed by the lines / ~ ~ ~ /C ¼ 0 and /f ¼ /C . Several high-resolution schemes based on different flux limiters have been proposed in literature and a comprehensive study of these schemes have been carried out by Waterson and Deconinck [19]. These authors have shown that the best convection schemes may be encapsulated into a single class of generalised piecewise linear schemes, referred to as GPL-j, whose flux limiter maybe defined as,





cðrÞHR ¼ max 0; min ð2 þ aÞr;

1þj 1j rþ ;M 2 2



where j defines the unbounded linear scheme (also called the base scheme). This flux limiter for the high resolution scheme may now be converted to its equivalent NV form using Eq. (3) easily. We remark that while the discussions on the FL and NV approaches have been limited to one-dimensional uniform meshes, the ideas may be extended to unstructured meshes (Fig. 1(b)) as well with relative ease. This can be achieved by realising that the value for the volume fraction in the far upwind cell (U) on unstructured meshes may be approximated as /U ¼ /D  2r/C  rCD [12]. 4. Design principles for interface capturing schemes While high-resolution schemes lead to reduced numerical diffusion as compared to their first-order counterparts, they still do not lead to sufficiently sharp resolution of discontinuities such as fluid–fluid interfaces that arise in multi-fluid flows. In order to maintain the interface sharp, it is therefore necessary to introduce further anti-diffusion in a controlled manner. The basic philosophy behind the design of robust schemes for interface capture is to combine high-resolution (HR) and compressive (BD) schemes using a suitable blending function [6]. There are however no guidelines provided on how these constituent schemes need to be chosen and their impact on the resulting interface capture scheme. In this section, we attempt to elaborate on this key philosophy, along the lines of the unified framework for higher-order schemes proposed in [19], to evolve a set of detailed guidelines for interface capturing schemes.

   r/  r  CD   f hf ¼ arccos   jr/f jjrCD j 2. Boundedness and linearity preservation: The interface capturing scheme must be bounded so that the values of volume fraction lie in [0, 1]. This means that the constituent HR and BD schemes must respect the CBC to ensure no undershoots and overshoots. Furthermore, in order to ensure second-order accuracy in regions of smooth solutions (r = 1), the HR scheme must be linearity preserving which implies that it must pass through (1/2, 3/4) in the NV diagram. 3. High-resolution scheme: The high-resolution scheme allows higher-order accuracy with monotonicity and is constructed using the following limiter, very similar to the GPL-j limiter in [19], where w refers to a constant quantity that may depend on the Courant number.





   1þj 1j ;M rþ 2 2

cðrÞHR ¼ max 0; min ð2 þ aÞr; w

4. Compressive scheme: The compressive scheme is responsible for sharpening the interface with its contribution being controlled by the blending function. There are no definitive strategies as yet to choose the BD scheme, and therefore we recommend that they be derived from the HR scheme itself by using the following limiter.

cðrÞBD ¼ max½0; minfð2 þ aÞr; Mg It is quite easy to see that the limiter has been derived from that of GPL-j scheme, by simply dropping the base unbounded linear scheme. 5. Influence and extent of downwinding: The concept of downwinding is related to anti-diffusion and determines the ability to sharpen the interface. The parameter a denotes the proximity to full downwinding while the maximum bound M determines the extent of downwinding from the compressive (BD) scheme. We choose a value of a ¼ 0 to ensure full downwinding, but its extent (value of M) will be determined based on numerical experiments. We remark that although different value of upper bounds for the HR and BD limiters maybe used, we recommend that the same value be used in both limiters. We remark here that while the blending function and properties of boundedness and linearity preservation for convection schemes are well-established, the construction of the limiters for high-resolution and compressive schemes, with a view to unify interface capturing schemes, is a novel proposal of the present study. Along with the choice of a suitable downwind parameter, these set of guidelines provide a complete unified framework that

J.K. Patel, G. Natarajan / Computers & Fluids 106 (2015) 108–118

can be used to analyse/develop new schemes for interface capture in multi-fluid flows.

111

dependence, defined by w ¼ ð1  CoÞ where Co is the local Courant number. A simple analysis also shows that the unbounded linear scheme corresponds to the QUICK scheme so that j ¼ 12. The HYPER-C is a very compressive scheme with the downwind extent parameter depending on Courant number. A close examination of the compressive limiter by using FL/NV transformations shows 1

that M = 2 Co  1 and a ¼ 0. The limiter for the BD scheme is clearly derived from that of the HR scheme. An important distinction of the CICSAM schemes from new schemes proposed in this work is the explicit dependence of the limiters on the Courant number. Clearly, for lower values of Courant number, a greater extent of downwinding is achieved along with a higher-order accurate convection scheme which confirms the superior performance of the scheme for Co < 0:3. For higher values of Courant number, the HR scheme becomes heavily biased to first-order upwind and the downwind influence diminishes considerably, which explains the loss of resolution and diffusion of the interfaces for higher Courant numbers. The performance of CICSAM may also be influenced, albeit to a lesser extent, by the choice of the blending function. The M-CICSAM scheme proposed in [24] is a definite improvement over the CICSAM scheme not just because it chooses a different HR scheme but also because it avoids Courant number dependence in the limiters employed for the compressive and high-resolution schemes. We refrain from the analysis of the M-CICSAM scheme, but remark that it shares nearly all characteristics of the FBICS scheme [12] which is discussed in the following, except for a different choice of the downwind influence parameter M and the blending function.

5. Investigations of existing interface capturing schemes This section is devoted to the analysis of existing interface capture schemes based on the unified framework derived using the design guidelines discussed in the previous section. In particular, we examine four different schemes viz. CICSAM, HRIC, InterGamma and FBICS schemes, starting from the NV formulation of their constituent HR and BD schemes. Using the FL/NV transformation, these schemes are cast as a class of GPL-j schemes, by deriving their specific limiters, thereby de-constructing how these schemes could have been possibly derived. The merits and limitations of these schemes are then discussed on the basis of the design guidelines which explain their success and potential failures. The FL and NV forms of these schemes and the specific flux limiters are summarised in Table 1. 5.1. Compressive Interface Capturing Scheme for Arbitrary Meshes, CICSAM [5] The CICSAM scheme was proposed by Ubbink and Issa for multi-fluid flows. The scheme switches (with a different blending function than discussed earlier) between the HYPER-C scheme and ULTIMATE QUICKEST scheme, the latter being the HR scheme and the former the BD scheme, and both satisfying CBC. The HR scheme preserves linearity as well, and employs for itself a blend of upwind and QUICK schemes with a Courant number Table 1 Schemes in flux limiters and normalised form. Schemes

Flux limiters

Values of M; w and

h n oi

1 1 ¼ max 0; min 2r; ð1  CoÞ 3r  1Þ 4 þ 4 ; 2ðCo

CICSAM [5]

cðrÞ cðrÞBD ¼ max½0; minf2r; 2ðCo1  1Þg

M ¼ 2ðCo1  1Þ w ¼ ð1  CoÞ j ¼ 1=2

HRIC [4]

cðrÞHR ¼ 0  0:7Co 0:7Co cðrÞBD ¼ min 2 0:70:3 r; 2 0:70:3

zfflfflfflfflfflffl}|fflfflfflfflfflffl{ 0:7  Co M¼2 0:7  0:3 w¼0

Inter-Gamma [9]

–

–

FBICS [12]

cðrÞHR ¼ max 0; min 2r; 2r þ 12 ; 4 cðrÞBD ¼ max½0; minf2r; 4g

M¼4 w¼1 j¼0

M-Gamma

cðrÞHR ¼ max½0; minfr; 3g cðrÞBD ¼ max½0; minf2r; 3g

M¼3 w¼1 j¼1

CUIBS

cðrÞHR ¼ max 0; min 2r; 23 r þ 13 ; 4 cðrÞBD ¼ max½0; minf2r; 4g

HR









K





M¼4 w¼1 j ¼ 1=3

j

Normalised form n ~ o ( ~ C þ3Þ / ~ BD ~C 6 1 min 8Co/C þð1CoÞð6 ;/ 06/ 8 f /~f HR ¼ ~C ~C P 1 / 0 P / n ~ o ( ~C  1 min 1; /CoC 0 2 > > >1 1 ~ C 6 1 and Co < 0:3 > 6 / > 2 > < 1 ~ ~ 0 < / 6 ð1 þ KÞ / C C 2 and 0:3 < Co < 0:7 /~f BD ¼ ~ C 6 1 and 0:3 < Co < 0:7 ~C 1 6 / > K þ ð1  KÞ/ > > 2 > > ~ C  1 and Co < 0:7 ~ >/ 0P/ > : C ~C ~ C and Co > 0:7 / for all values of / 8 ~ 2 þ 3/ ~C 0 < / ~C 6 1 > < 2/ c 2 1 ~ /~f BD ¼ 1 2 < /C 6 1 > :~ ~ 0 P /C P 1 /C 8 ~C 6 1 ~C > 0 < / 3 / > 8 > <~ ~C 6 3 /C þ 14 18 < / 4 /~f HR ¼ 3 ~ > 1 > 4 < /C 6 1 > :~ ~C P 1 0P/ 8 /C ~ ~ C  1=3 < 3/C 0 < / ~C  1 /~f BD ¼ 1 1=3 < / :~ ~C P 1 /C 0P/ 8 5~ ~C 6 1 > 0 :~ ~ 0 P /C P 1 8 /C ~C 6 2 ~C 0 < / > 52 / < 5 2 ~ /~f BD ¼ 1 5 < /C  1 > :~ ~ 0 P /C P 1 /C 8 ~C 6 2 ~C > 0 13 >5 < ~C þ 1 2 < / ~C 6 4 / 3 13 5 /~f HR ¼ 6 4 ~ > 1 < / 6 1 > C 5 > :~ ~ 0 P /C P 1 8 /C ~C 0 < / ~C 6 1 > 3/ < 3 BD ~ 1 ~C  1 /f ¼ 1 < / 3 > :~ ~C P 1 /C 0P/

112

J.K. Patel, G. Natarajan / Computers & Fluids 106 (2015) 108–118

5.2. High Resolution Interface Capturing scheme, HRIC [4] The HRIC scheme attempts to avoid explicit Courant number dependence that plagues CICSAM scheme and to simplify the calculation of face values for the variables, which in the present study is the volume fraction. The overall scheme is characterised by a three-step method to enforce a switching in Courant-number domain. Unfortunately, this brings in dependence of the Courant number into the HRIC scheme, which in fact it sought to eliminate. Moreover, the use of FL/NV transformations to each of the three steps indicates that the HR scheme for very low (less than 0.3) and very high (greater than 0.7) Courant numbers is merely the first-order upwind scheme, with a good amount of compression in the former case which deteriorates to zero in the latter case. It therefore follows that w ¼ 0 and the value of M = 2 (for Co < 0:3) and M = 0 (for Co > 0:7). For intermediate values of Courant number 0:3 < Co < 0:7, the HR scheme remains first-order upwind with w ¼ 0, while the compressive scheme involves downwinding whose extent and influence are dependent on Courant number, reducing as Courant number increases. Specifically, a closer examination using FL/NV formulations show that a ¼ 2K  2 and M = 2K 0:7Co where K ¼ 0:70:3 is a Courant-number dependent constant. For all three Courant number regimes, the lack of any genuinely high-order accurate scheme in the blending makes HRIC a more diffusive scheme than CICSAM, and the compressive scheme only serves to selectively diminish this over-diffusion. We remark that the scheme is fully first-order upwind for high Courant numbers and remains robust at the expense of severe loss of resolution of interface sharpness. The fact that the HR scheme for HRIC is infact a low-resolution scheme explains the superior performance of CICSAM in comparison to HRIC in the studies in [25]. The HRIC scheme in its original form would be advisable only on very fine meshes with a severe Courant number restriction if satisfactory results are to be obtained. 5.3. Inter-Gamma scheme [9] The Inter-Gamma was proposed to specifically target sharp resolution of discontinuities by introducing significant downwinding within a single advection scheme. Consequently, it involves no blend of HR and BD schemes but rather behaves merely as a BD scheme. There is no control on amount and manner of downwinding that must be applied in response to the orientation of the interface and the NV formulation of the scheme (Table 1) shows a quadratic segment which violates linearity preservation. The poor performance of the Inter-Gamma scheme for simple advection problems in [8] may be attributed to these characteristics of the scheme, although it is shown to work well for some practical flow problems with a severe Courant number restriction. The authors are of the opinion that any scheme that violates the requirements outlined in Section 4 does not qualify as a ‘‘good’’ interface capturing scheme and its excellent performance for some practical problems is rather fortuitous. 5.4. Flux Blending Interface Capturing Scheme, FBICS [12] Tsui and co-workers have developed two novel schemes for interface captures based on flux blending, FBICS-A and FBICS-B and compared their performance with HRIC and CICSAM schemes for typical advection problems. The FBICS-A scheme (hereafter referred to simply as FBICS) was found to be most efficient in capturing the discontinuities with low errors that were not dependent on the Courant number. FBICS is built to satisfy CBC, and uses the Fromm scheme as the base linear scheme. The study in [12] however does not discuss the construction of the limiters employed for the HR and BD schemes.

It is easy to see that the HR scheme is linearity preserving and the FL/NV transformations show that the choice of j ¼ 0 and M = 4 would result in the flux limiter for the HR scheme. The limiter for BD scheme may then be easily derived by dropping the unbounded Fromm scheme. The specific choice of the value of M remains unexplained in the original work and we shall attempt an explanation of this choice in the subsequent sections. It may be remarked that if M = 2, the FBICS scheme would reduce to the M-CICSAM scheme, with the only difference then being the choice of blending function whose effect is less significant in the overall performance. The authors believe that while the schemes have been independently developed, the present analysis show that the two schemes are only minor variants of each other. 6. Development of new schemes for interface capture The design principles and unified framework are now employed to develop robust and accurate interface capturing schemes devoid of complex geometric reconstruction. Two different schemes are constructed, viz. the M-Gamma scheme and CUIBS, and the former may be considered as an improved version of the Inter-Gamma scheme presented earlier. These schemes are designed to work with arbitrary polygonal meshes and as will be demonstrated in later sections, exhibit a performance that is nearly independent of the Courant number unlike HRIC and CICSAM schemes. The FL/ NV description of the limiters used to construct the HR and BD schemes for these novel interface capture schemes are summarised in Table 1. 6.1. Modified Gamma scheme (M-Gamma) The failure of the Inter-Gamma scheme motivates the construction of a new interface capturing scheme within the class of blended GPL-j schemes that respects the design guidelines. We begin by considering the original Gamma scheme [23], which is a HR scheme that in itself combines central differencing scheme with an ad hoc quadratic variation, controlled by a parameter b. ~ C 2 ½0; b is The quadratic variation which is applied in the range / defined by,

 ~2  ~C ~ f ¼  /C þ 1 þ 1 / / 2b 2b Comparing this variation with the quadratic variation adopted in Inter-Gamma scheme (Table 1), it follows that b ¼ 1=4. We construct the HR scheme for M-Gamma simply by considering a linear ~ C 2 ½0; 1=4 followed by the central differencvariation in the range / ing scheme (which is the base linear scheme). The HR scheme may be considered as a special case of GPL-j schemes where j ¼ 1 and M = 3 with the BD scheme obtained by discarding the base scheme, resulting in the limiters shown in Table 1. The value of M is fixed here by the choice of b as can be seen from a closer consideration of the NV diagram shown in Fig. 2(a) which also indicates that the HR scheme is linearity preserving. A different choice for M in both HR and BD limiters may be employed (such as that of FBICS where M = 4) or the values of M in the limiters can be differentially selected but these ideas are not explored in this study. 6.2. Cubic Upwind Interpolation based Blending Scheme (CUIBS) The CUIBS scheme is directly inspired by the work of Waterson and Deconinck [19] which shows that the best convective schemes belong to the GPL-j class of schemes. These schemes are MUSCL, SMART and Koren’s limited CUI scheme and perform excellently on smooth and discontinuous advection problems. While the first two schemes have been employed as HR schemes for interface

113

J.K. Patel, G. Natarajan / Computers & Fluids 106 (2015) 108–118

Fig. 2. Normalised variable diagram (a) M-Gamma and (b) CUIBS.

capture to varying extent in the past (in FBICS and STACS [6] for instance), there have been no attempts to date to use the limited CUI scheme as a building block for interface capturing schemes. We merely use the Koren limited CUI scheme as the HR scheme, which corresponds to a HR limiter with j ¼ 1=3 and M = 4. We remark however that the M value is different from the original CUI scheme and the choice of this value is justified through numerical experiments in the following section. The BD limiter for the compressive scheme may be obtained by dropping the base scheme and the resulting limiter is same as that used for FBICS scheme. The HR scheme is linearity preserving and the normalised variable diagram for this scheme is depicted in Fig. 2(b). It is expected that the superior characteristics of this unique thirdorder upwind scheme in the HR limiter will result in CUIBS having improved interface capture performance in comparison to its counterparts. We emphasise that this work is the first instance of the CUI scheme being employed for sharp resolution of interfaces. 7. Investigations on advection problems The new schemes proposed in this study are tested for the ability to capture interfaces and preserve their sharpness using three advection problems. Tests are carried out by solving the advection equation alone with prescribed velocity fields. The M-Gamma and CUIBS schemes are also compared with the FBICS scheme, to understand the relative merits of different schemes. A computational domain of 4  4 is chosen and the Courant number and error are defined as,

P f 2c maxðuf  Dsf ; 0ÞDt Co ¼ Xc Pnc n e c¼1 k/c Xc  /c Xc k E¼ Pnc 0 c¼1 /c Xc where nc is the number of volumes in the domain, c denotes the cell, f denotes the faces of the cell, Dsf represents the vectorial area and superscripts n; e and 0 represent the numerical and exact solutions at the final time level and the exact solution at the initial time level respectively. Tests are performed for small, intermediate and large Courant numbers on structured and unstructured meshes to extensively analyse the proposed interface capturing schemes.

mesh and 22,584 unstructured mesh. Table 2 shows the contour plots of hallow square after 1 unit time on structured and unstructured mesh employing FBICS, M-Gamma and CUIBS schemes. Simulation is executed for low and high Courant number (Co) 0.15, and 0.8 respectively aimed to check the Courant number dependency of present implicit methodology. The performance of FBICS is in line with the work of Tsui et al. [12] wherein similar results have been reported. Despite of slight diffusive behaviour on structured mesh, M-Gamma scheme produces better results on unstructured mesh by virtue of finer mesh. Performance of CUIBS scheme is as good as FBICS scheme on both mesh arrangements. Error tendency with Courant number are reported in Fig. 3(a) and (b). 7.2. Test 2: Advection of a circle in shear flow The test consists of a circle of radius 0.2p units centred at (1.5, 0.8) and subjected to a velocity field ðu; v Þ = (sin(x)cos(y), cos(x)sin(y)). The circle is strained for N time units after which the velocity is reversed and the circle returns to its original configuration after further N time units. Computations are done on 100  100 structured mesh and a triangular mesh consisting of 22,584 cells. The contours of the volume fraction for N = 8 for the three schemes on both meshes for Co = 0.8 are shown in Table 3. The error variation with Courant number for these grids are depicted in Fig. 4. While all three schemes show Courant number independent errors, the M-Gamma scheme is found to be inferior to the FBICS and CUIBS schemes. However, M-Gamma is no longer restricted by a Courant number restriction unlike the Inter-Gamma scheme. While the FBICS and CUIBS schemes are comparable on structured meshes, the latter outperforms on unstructured grid in terms of the error. Table 2 /-contours for Test 1. Contour levels are 0.05:0.05:0.95. Schemes

FBICS [12]

M-Gamma

7.1. Test 1: Advection of hollow square in oblique flow To validate the accuracy of interface capturing schemes, hollow square of 0.8 outer width and 0.4 inner width and centred at ð0:8; 0:8Þ, is advected in uniform oblique velocity field of ðu; v Þ ¼ ð2; 1Þ. Simulation is carried out for 100  100 structured

CUIBS

Structured mesh

Unstructured mesh

Co = 0.15

Co = 0.15

Co = 0.8

Co = 0.8

114

J.K. Patel, G. Natarajan / Computers & Fluids 106 (2015) 108–118

(b) 0.28

FBICS M-Gamma CUIBS

Numerical error

0.26 0.24 0.22 0.2 0.18 0.16

0.28

Numerical error

(a)

FBICS M-Gamma CUIBS

0.26

0.24

0.14 0.12 0.1 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.1

0.2

Courant number (Co)

0.3

0.4

0.5

0.6

0.7

0.8

Courant number (Co)

Fig. 3. Numerical error against Courant number for hollow square using (a) structured mesh and (b) unstructured mesh.

Table 3 /-contours for Test 2 at Co = 0.8. Contour levels are 0.05:0.05:0.95. Schemes

Structured mesh Forward

Table 4 /-contours for Test 3. Contour levels are 0.05:0.05:0.95.

Unstructured mesh Backward

Forward

Schemes

Backward

FBICS [12]

Structured mesh

Unstructured mesh

Co = 0.15

Co = 0.15

Co = 0.8

Co = 0.8

FBICS [12]

M-Gamma M-Gamma

CUIBS CUIBS

7.3. Test 3: Zalesak solid body rotation This test consists of the advection of a slotted disk in a rotational flowfield. The disk has unit diameter and is centred at (2, 2.75) and is cut by a slot of 0.12 unit width. The disk is rotated around centre of rotation (x0 ; y0 Þ ¼ ð2; 2Þ with a constant angular velocity of 0.5. The velocity field is therefore ðu; v Þ ¼ ðxðy  y0 Þ; xðx  x0 ÞÞ. Simulations are carried out for one rotation using structured mesh of 200  200 and unstructured mesh of 40,532

(a)

(b)

0.14 0.13

0.1 0.09 0.08 0.07

0.11 0.1 0.09 0.08 0.07

0.06

0.06

0.05

0.05

0.04 0.1

FBICS M-Gamma CUIBS

0.12

Numerical error

0.11

0.14 0.13

FBICS M-Gamma CUIBS

0.12

Numerical error

triangles. Table 4 shows the contours of volume fraction for different Courant numbers and the error variation on both grids as a function of Courant number is shown in Fig. 5. While this case gives further evidence of performance not depending on Courant number, the M-Gamma and CUIBS schemes are found to have lower errors compared to the FBICS scheme on unstructured meshes. On structured meshes, all three schemes show very similar behaviour. The contour plots indicate that the FBICS scheme distorts the profiles on unstructured meshes as compared to the new schemes proposed here. This test case was not investigated in [12] and this observation is surprising. The CUIBS scheme is at least as competitive as FBICS, and the improved performance of

0.04 0.2

0.3

0.4

0.5

0.6

Courant number (Co)

0.7

0.8

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Courant number (Co)

Fig. 4. Variation of numerical error with Courant number for Test 1 (a) structured mesh and (b) unstructured mesh.

0.8

115

J.K. Patel, G. Natarajan / Computers & Fluids 106 (2015) 108–118

(a)

(b)

0.4 FBICS M-Gamma CUIBS

0.3 0.25 0.2 0.15 0.1 0.05

FBICS M-Gamma CUIBS

0.35

Numerical error

Numerical error

0.35

0.4

0.3 0.25 0.2 0.15 0.1 0.05

0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0 0.1

0.2

0.3

Courant number (Co)

0.4

0.5

0.6

0.7

0.8

Courant number (Co)

Fig. 5. Variation of numerical error with Courant number for Test 2 (a) structured mesh and (b) unstructured mesh.

0.8

−4

M M M M M

0.7

0.5 0.4 0.3 0.2 0.1 0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

x 10

3.085

Mass error (Emass)

Numerical error

0.6

3.09

=2 =3 =4 =6 = 10

3.08 3.075 3.07 3.065 3.06

0.8

Courant number (Co) 3.055 Fig. 6. Error variation with Courant number for different values of the upper bound M on unstructured mesh.

the former for some cases maybe attributed to the third-order accurate base scheme and the marginally lesser extent of downwinding introduced in the constituent HR scheme. Nevertheless, the behaviour of M-Gamma scheme is somewhat perplexing, and it is conjectured that the downwinding extent plays a major role. In particular, the upper bound M may have an optimal value that depends on the test problem and/or grid employed which decides the superior interface capturing scheme for a particular problem and grid topology. We precisely investigate the role of the maximum bound, albeit in a limited sense, in the following numerical experiment. 7.4. Effect of extent of downwinding We study the effect of the maximum bound M in the flux limiters of the schemes to understand its role in interface capture and determine an optimal value, if possible. We consider the advection of a hollow square translated obliquely by a constant velocity field ðu; v Þ ¼ ð2; 1Þ on unstructured grids using the CUIBS scheme. The error variation with Courant number plotted in Fig. 6 for different values of M show that the errors are minimal for M 2 [2,4]. This suggests the existence of a possibly ‘‘optimal’’ range for M which may be understood from the fact that smaller M results in lesser downwinding from the compressive scheme and hence more diffusion while a higher value of M may result in over-compression.

0

4

8

12

16

Non−dimensional time (t) Fig. 7. Relative change of the total mass of the volume fraction in the shearing flow test using CUIBS scheme on unstructured mesh.

This study justifies M = 4 as a good choice in terms of error and explains its choice in CUIBS and FBICS schemes. We remark, as also discussed earlier, that a different value of M could become optimal for the schemes in certain interface problems and grid topologies and therefore the value of M = 4 is by no means a universal choice. 7.5. Mass conservation An important and often overlooked issue in interface capturing schemes and VOF approaches, particularly on unstructured meshes, is the conservation of mass of the fluids. Despite solving the conservation laws in a conservative framework, the clipping of the volume fraction (to ensure it is bounded between 0 and 1 in each cell) may introduce errors that can violate discrete conservation. Every fluid in the domain has a constant density and their mass in a closed domain must be conserved over time and therefore we quantify the errors in discrete mass conservation as,

Pnc Emass ¼

n 0 c¼1 /c c  /c Pnc 0 c¼1 /c c

X

X

Xc

ð4Þ

where the subscript and superscript denote the volume and time level where the volume fraction is computed.

116

J.K. Patel, G. Natarajan / Computers & Fluids 106 (2015) 108–118

(b)

0.4

Displacement from mean position

Displacement from mean position

(a)

Present−CUIBS −−− Tezduyar et. al. [26]

0.3 0.2 0.1 0 −0.1 −0.2

0

10

20

30

40

0.3 Present−CUIBS −−− Tezduyar et. al. [26]

0.2 0.1 0 −0.1 −0.2 −0.3 −0.4

Non−dimensional time (t)

0

10

20

30

40

Non−dimensional time (t)

Fig. 8. The time history of interface location relative to the average height of 0.3 on the (a) left wall and (b) right wall.

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.2

0.4

0.6

0.8

0

0

0.2

(a) t=0 0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.2

0.4

0.4

0.6

0.8

0.6

0.8

(b) t=5

0.6

0.8

(c) t=20

0

0

0.2

0.4

(d) t=40

Fig. 9. Two fluid interface using CUIBS on structured mesh along with velocity vector.

To investigate the mass loss (or gain) due to limiting of the volume fraction, especially on non-orthogonal meshes with time, we consider the advection of circle in shear flow as discussed in Section 7.2 using the CUIBS scheme on unstructured meshes. The mass error Emass is evaluated for this case and its temporal history is shown in Fig. 7. We see that the average mass error is around 104 and remains bounded, ensuring that the present approach satisfies the conservation criterion to engineering accuracy. The loss/gain in mass is therefore negligible and is not expected to significantly affect the solutions to complex interfacial flows, as is demonstrated in the test cases considered in Section 8.

8. Application to binary immiscible fluid flows The proposed CUIBS scheme, which has shown the best performance among all schemes discussed herein, is applied to engineering problems of interest involving binary interfacial flows to assess its utility. The governing Navier–Stokes equations are solved to obtain the velocity field using a hybrid staggered/non-staggered approach with a fractional step method, that is second-order accurate in space and time. The pressure field is determined by solving a variable-coefficient Poisson equation in an implicit fashion, akin to the solution of the volume fraction equation. A single fluid for-

117

J.K. Patel, G. Natarajan / Computers & Fluids 106 (2015) 108–118

4

4

4

4

4

3.5

3.5

3.5

3.5

3.5

3

3

3

3

3

2.5

2.5

2.5

2.5

2.5

2

2

2

2

2

1.5

1.5

1.5

1.5

1.5

1

1

1

1

1

0.5

0.5

0.5

0.5

0.5

0 -0.5

0

0.5

0 -0.5

0

0.5

0 -0.5

0

0.5

0 -0.5

0

0.5

0 -0.5

0

0.5

Fig. 10. Rayleigh–Taylor instability: contours of volume fraction = 0.5 using CUIBS scheme at times 0.5, 1, 1.5, 2, 2.5.

8.1. Two fluid interface problem We consider the viscous sloshing in a 0.8  0.6 domain with two fluids, with the lighter fluid above and heavier below and the initial interface defined by y = 0.4–0.25x. No-slip condition is enforced on the horizontal top and bottom walls while the vertical walls are defined by a velocity slip condition. The Reynolds number pffiffiffiffiffi and Froude number defined by Re ¼ gLLqmin =lmin and Fr ¼ U 2 =gL are chosen to be 300 and 0.294 respectively, where the lighter fluid is taken as the reference. The density ratio is taken as 2 and the viscosities of the two fluids are considered equal. The domain is divided into 24,220 rectangular elements and the results of the simulations are compared with those obtained by Tezduyar et al. [26]. Fig. 8 shows the temporal history of the interface location on the left and right walls, which are in good agreement with the FEM results in [26] and show a damped behaviour indicating a steady state at very long times. Fig. 9 shows the interface at four different time instants (t = 0, 5, 20, 40) indicating that the initial sharp resolution of the interface is preserved to a great extent by the CUIBS scheme. 8.2. Rayleigh–Taylor instability The R–T instability, which arises when a heavier fluid is on top of a lighter one, poses a more challenging problem for interface capturing schemes owing to the complex dynamics of the interfacial instability. We consider a 1  4 domain is discretised into 21,456 triangular elements, with heavier liquid on the top and lighter liquid on the bottom such that the viscosities are equal and the density ratio is 3. This corresponds to an Atwood number qmin Þ At ¼ ððqqmax  ¼ 0:5 and the initial perturbed interface is defined by max þqmin Þ y ¼ 2 þ a cosð0:5pxÞ with amplitude of a = 0.05. This problem is known to be quite grid sensitive and therefore an unstructured grid symmetric about the vertical centreline is used to minimise

2.8

Displacement from mean position

mulation is employed wherein the face values of density and viscosity are obtained from the centroidal values of these scalars using simple and harmonic averaging respectively. In the studies discussed below, only binary immiscible fluids with low density and viscosity ratios (typically lesser than 10) are considered although the extensions to high density ratios and multiple fluids are relatively straightforward.

Present−CUIBS −−− Guermond et. al. [28] ooo Tryggvason [27]

2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8

0

0.5

1

1.5

2

2.5

Non−dimensional time = t √(At) Fig. 11. Position of rising and falling bubbles versus time on unstructured meshes using CUIBS scheme.

any asymmetry errors. Effects of surface tension are neglected and the Reynolds and Froude numbers are set to 1000 and 1 respectively. Fig. 10 p reveals the time history of the interface at different time ffiffiffiffiffi levels of t At ¼ 0:5; 1; 1:5; 2; 2:5. The position of the rising bubble at the side wall and falling bubble at the mid section of are plotted against the non-dimensional time in Fig. 11, and are compared with numerical simulations of Trggvasaon [27] and Guermond and Quartapelle [28]. The agreement of the results from the present simulations with those in [27,28] is excellent, and further validates the robustness and efficacy of the CUIBS scheme to handle complex flow phenomena in realistic fluid flows.

9. Conclusions A generic and unified framework for analysis and construction of interface capturing schemes has been developed in this work. The framework exploits the key philosophy of blending compres-

118

J.K. Patel, G. Natarajan / Computers & Fluids 106 (2015) 108–118

sive and high-resolution schemes using a class of generalised piecewise limiters and a definitive set of design principles for their construction is proposed. This framework is used to analyse four existing schemes and their success and/or failure are discussed based on the set of design principles outlined. Two novel schemes, M-Gamma and CUIBS are proposed on the basis of this framework and are shown to be bounded and stable through studies on typical advection problems. These studies also demonstrate that the proposed schemes track the interface sharply on structured and unstructured grids, and that their performance has no significant dependence on the Courant number. Numerical investigations of viscous sloshing and Rayleigh–Taylor instability using the CUIBS scheme, which shows the best performance among all schemes for typical advection tests, are performed to assess the scheme on different grid topologies. While the studies herein have been limited to density-stratified flows, the guidelines and unified framework presented in this study are generic and apply equally well to viscosity-stratified flows as well. The CUIBS scheme is found to preserve sharp interfaces without distortion while resolving interfacial phenomena with good accuracy, making it a promising choice for interface capture in density and viscosity stratified multiphase flows on three-dimensional unstructured meshes. References [1] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39:201–25. [2] Noh WF, Woodward P. SLIC (simple line interface calculations). Lect Notes Phys 1976;59:330–40. [3] Youngs DL. Time-dependent multi-material flow with large fluid distortion. In: Morton K, Baines M, editors. Numerical methods for fluid dynamics. New York: Academic Press; 1982. p. 273–85. [4] Muzaferija S, Peric M, Sames P, Schellin T. A two-fluid Navier–Stokes solver to simulate water entry. In: Proceeding of twenty-second symposium on naval hydrodynamics, Washington, DC; 1998. p. 638–49. [5] Ubbink O, Issa RI. A method for capturing sharp fluid interfaces on arbitrary meshes. J Comput Phys 1999;153:26–50. [6] Darwish M, Moukalled F. Convective schemes for capturing interfaces of freesurface flows on unstructured grids. Numer Heat Transfer, Part B 2006;49:19–42. [7] Moukalled F, Darwish M. Transient schemes for capturing interfaces of freesurface flows. Numer Heat Transfer, Part B 2012;61:171–203. [8] Gopala VR, van Wachem BGM. Volume of fluid methods for immiscible-fluid and free-surface flows. Chem Eng J 2008;141:204–21.

[9] Jasak H, Weller HG. Interface-tracking capabilities of the Inter-Gamma differencing scheme. Technical report, Imperial College, University of London; 1995. [10] Walters DK, Wolgemuth NM. A new interface-capturing discretization scheme for numerical solution of the volume fraction equation in two-phase flows. Int J Numer Methods Fluids 2009;60:893–918. [11] Xiao F, Honma Y, Kono T. A simple algebraic interface capturing scheme using hyperbolic tangent function. Int J Numer Methods Fluids 2005;48:1023–40. [12] Tsui YY, Lin SW, Cheng TT, Wu TC. Flux-blending schemes for interface capture in two-fluid flows. Int J Heat Mass Transfer 2009;52:5547–56. [13] Heyns JA, Malan AG, Harms TM, Oxtoby OF. Development of a compressive surface capturing formulation for modelling free-surface flow by using the volume-of-fluid approach. Int J Numer Methods Fluids 2013;71:788–804. [14] Osher S, Sethian LA. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79:12–49. [15] Chessa J, Belytschko T. An enriched finite element method and level sets for axisymmetric two-phase flow with surface tension. Int J Numer Methods Eng 2003;13:2041–64. [16] Sahu KC, Vanka SP. A multiphase lattice Boltzmann simulations of buoyancy induced mixing in a tilted channel. Comput Fluids 2011;50:199–215. [17] Redapangu PR, Vanka SP, Sahu KC. A multiphase lattice Boltzmann simulations of buoyancy-induced flow of two immiscible fluids with different viscosities. Eur J Mech – B/Fluids 2012;34:105–14. [18] Kim J. A diffuse-interface model for axisymmetric immiscible two-phase flow. Appl Math Comput 2005;160:589–606. [19] Waterson NP, Deconinck H. Design principles for bounded higher-order convection schemes – a unified approach. J Comput Phys 2007;224:182–207. [20] http://www.ssisc.org/lis; 2013. [21] Roe PL. Characteristic-based schemes for the Euler equations. Annu Rev Fluid Mech 1986;18:337. [22] Gaskell PH, Lau AKC. Curvature-compensated convective transport: SMART, a new boundedness-preserving transport algorithm. Int J Numer Methods Fluids 1988;8:617–41. [23] Jasak H, Weller HG, Gosman AD. High resolution NVD differencing scheme for arbitrarily unstructured meshes. Int J Numer Methods Fluids 1999;31:431–49. [24] Waclawczyk T, Koronowicz T. Remarks on prediction of wave drag using VOF method with interface capturing approach. Arch Civ Mech Eng 2008:8. [25] Waclawczyk T, Koronowicz T. Comparison of CICSAM and HRIC highresolution schemes for interface capturing. J Theor Appl Mech 2008;46:325–45. [26] Tezduyar TE, Behr M, Mittal S. A new strategy for finite element computations involving moving boundaries and interfaces – the deforming-spatial-domain/ space-time procedure: II. Computation of free-surface flows, two-liquid flows, and flows with drifting cylinders. Comput Methods Appl Mech Eng 1992;94:353–71. [27] Tryggvason G. Numerical simulation of the Rayleigh–Taylor instability. J Comput Phys 1988;75:253–82. [28] Guermond JL, Quartapelle L. A projection FEM for variable density incompressible flows. J Comput Phys 2000;165:167–88.