The MAC method

Available online at www.sciencedirect.com

Computers & Fluids 37 (2008) 907–930 www.elsevier.com/locate/compfluid

Review

The MAC method S. McKee a,*, M.F. Tome´ b, V.G. Ferreira b, J.A. Cuminato b, A. Castelo b, F.S. Sousa b, N. Mangiavacchi c a

Department of Mathematics, University of Strathclyde, Livingstone Tower, 26 Richmond Street, Glasgow G1 1XH, UK b Department of Applied Mathematics and Statistics, University of Sa˜o Paulo, USP, Sa˜o Carlos, SP, Brazil c Department of Mechanical Engineering, State University of Rio de Janeiro, UERJ, Rio de Janeiro, Brazil Received 1 January 2007; accepted 4 October 2007 Available online 24 November 2007

Abstract In this article recent advances in the Marker and Cell (MAC) method will be reviewed. The MAC technique dates back to the early 1960s at the Los Alamos Laboratories and this article starts with a historical review, and then a brief discussion of related techniques. Improvements since the early days of MAC (and the Simplified MAC – SMAC) include automatic time-stepping, the use of the conjugate gradient method to solve the Poisson equation for the corrected velocity potential, greater efficiency through stripping out the virtual particles (markers) other than those near the free surface, and more accurate approximations of the free surface boundary conditions, the addition of bounded high accuracy upwinding for the convected terms (thereby being able to solve higher Reynolds number flows), and a (dynamics) flow visualization facility. More recently, effective techniques for surface and interfacial flows and, in particular, for accurately tracking the associated surface(s)/interface(s) including moving contact angles have been developed. This article will concentrate principally on a three-dimensional version of the SMAC method. It will eschew both code verification and model validation; instead it will emphasize the applications that the MAC method can solve, from multiphase flows to rheology. Ó 2007 Elsevier Ltd. All rights reserved.

Contents 1. 2. 3. 4. 5.

6.

*

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Background of the marker and cell method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A review of interfacial tracking techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The emphasis of the review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Governing equations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1. Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2. Interfacial and free surface conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Discretization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1. Markers and cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2. Initial approximation of the free surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3. Interfacial tension effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4. Surface merging and splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4.1. Merging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4.2. Demerging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Corresponding author. Tel.: +44 0 141 548 3671; fax: +44 0 141 548 3345. E-mail address: [email protected] (S. McKee).

0045-7930/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2007.10.006

908 908 910 911 911 912 912 912 914 914 914 915 915 916

908

7.

8.

9.

10. 11.

12.

13.

14.

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

Numerical procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1. Time-stepping procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2. Finite difference approximation of the convective terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Undulation removal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.1. Normal vector calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2. Vertex balance procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3. Undulation removal procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Curvature calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.1. Normal vector approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.2. Curvature approximation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A dynamic contact angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Data structure and system environment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1. Data structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2. The module Modflow-3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.3. The module Simflow-3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.4. The module Visflow-3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.5. Communication between the modules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.1. Jet buckling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.2. Splashing drop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.3. Hydraulic jump . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.4. Rising bubbles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.5. Container filling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Viscoelastic flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13.1. Oldroyd-B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13.2. Elastic jet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13.3. Die-Swell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Future developments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

916 917 917 917 918 918 919 919 920 920 920 921 921 922 922 923 923 923 923 923 924 924 925 925 925 925 927 927 928 928

1. Introduction

2. Background of the marker and cell method

The purpose of this article is to provide an overview of the recent advances in the MAC method. Although it is a technology dating back to the early sixties at the Los Alamos Laboratories, with the considerably greater computing power that we now enjoy, it is witnessing a revival and, for free surface fluid flow problems, showing itself the equal of any of the competing methods. This review article will provide a brief history of the MAC technique, and even briefer description of related methods but will include a fairly complete description of the techniques and ideas behind the method. It will show that, with high order monotone upwinding approximation of the inertial terms, high Reynolds number flows are attainable and this will be illustrated by a simulation of the hydraulic jump. Viscous jet flow will be studied, particularly viscous jet buckling, examples of container filling will be given, and a liquid drop splashing onto a quiescent fluid will be simulated. Finally, we shall show how to extend the MAC method to viscoelastic fluid flow, thus opening up the field of computational free surface rheology, to which the method is particularly well adapted.

There was a prolific development of computational fluid dynamics (CFD) methods in the fluid dynamics group at Los Alamos Laboratories in the years from 1958 to the late 1960s. This development was largely due to the energy, creativity and leadership of Francis Harlow. It was he who in 1958 developed the original particle-in-cell (PIC) code [21,35], which used mass particles that carried material position, mass and species information. Despite the special capabilities of PIC for following discontinuities, it did not give accurate solutions in general because the transfer of information between the particles and the underlying grid resulted in numerical diffusion. Twenty years ago, in the 1980s, there was a resurgence of interest in particle methods, resulting in a variant called FLIP (e.g. see [7,8]); FLIP overcame this difficulty of numerical diffusion to a large extent. The principal applications of FLIP have been in magnetohydrodynamic flows, in the earth’s magnetosphere and the sun’s heliosphere (see e.g. [9]). The MAC method first appeared in 1965. It was developed by Harlow and Welch [36,104] specifically for free surface flows as a variant of PIC method. Based on an

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

Eulerian staggered grid system, the MAC method is a finite difference solution technique for investigating the dynamics of an incompressible viscous fluid. It employs the primitive variables of pressure and velocity. One of the key features is the use of Lagrangian virtual particles, whose coordinates are stored, and which move from a cell to the next according to the latest computed velocity field. If a cell contains a particle it is deemed to contain fluid, thus providing flow visualization of the free surface. In 1970, Amsden and Harlow [1] subsequently developed a simplified MAC method (SMAC) which circumvented difficulties with the original method by splitting the calculation cycle into two parts, namely: a provisional velocity field calculation followed by a velocity revision employing an auxiliary potential function to ensure incompressibility throughout. In this report, Amsden and Harlow describe a specific program, ZUNI, which embodies the SMAC methodology. This code was used to calculate two-dimensional flows in rectangular or cylindrical coordinates. It could deal with free surface flows and with free-slip or no-slip conditions applied on rigid boundaries. In addition, provision was made for prescribed inflows and outflows, and for the placing of a rectangular obstacle within the flow field. Fundamentally though, the fluid domain was required to be rectangular and it was restricted to two dimensions. The MAC method was the basis of the later ‘‘particleless” techniques for compressible and incompressible flows: many variants of the SMAC code appeared over the intervening years and are actively continuing to appear today. Due to the limitation on the time step-size, Deville [18], Pracht [72] and Golafshani [32] developed an implicit scheme. In particular, Pracht [72] developed an implicit treatment of the velocity similar to the implicit-fluid Eulerian (ICE) method (see [39]), using an arbitrary Lagrangian–Eulerian (ALE) computational mesh. This, in principle, allowed the calculation of flows involving curved or moving boundaries. For recent work employing these ideas, see [46]. Hirt and Shannon [41] and Nichols and Hirt [64] were concerned with a loss of accuracy on the free surface and suggested an improved treatment for the free surface stress conditions whereby the normal stress was applied on the actual fluid surface rather than on the centre of the surface cells. The first ideas to be distributed internationally were the SOLA family of codes written by Hirt, Nichols and others in the early 1970s. These include extensions to two immiscible fluids in the SOLA-VOF code [42,43]. One member of the SOLA family was the SOLA-DF code which included multiphase treatment (see [99]). About the same time, Butler and Rivard and others [10,75] developed the first reactive code, RICE; this evolved into the most widely used code, APACHE-CONCHAS-KIVA [74]. It is perhaps less well known that back in the 1960s Harlow and his co-workers also contributed to early numerical modeling of turbulence and, in particular, by the postulation of the j–e model, now well known world-wide [37,38]. Work in this area continues: see, for example, a

909

recent application of the SMAC method to the j–e equations by Ferreira et al. [24]. The 1980s saw many applications of the SMAC method: Sicilian and Hirt [79] developed a Containment Atmosphere Prediction (CAP) code using the PIC approach to model the flow of a containment; Miyata and Nishimura [61] and Miyata [62] used SMAC for the simulation of both water waves generated by ships and breaking waves over circular and elliptical bodies. At the same time, the UK electricity industry (then the Central Electricity Generating Board) was taking an interest in the MAC method. McQueen and Rutter [60] and Markham and Proctor [59] described modifications to the original fluid flow code ZUNI to provide enhanced performance. Essentially, they employed an automatic time-stepping routine and a preconditioned conjugate gradient solver for the Poisson equation in a rectangular region. In the 1990s many authors considered different, but related methods, like volume of fluid (VOF) [43], Euler–Lagrangian tracking methods and, more recently, level set methods, a survey about which is given in the next section. Nevertheless, there was a revival of the MAC method, principally by a research group in Brazil. Motivated initially by a research contract from Unilever Research plc to simulate container filling, Tome´ and McKee [88] substantially extended the work of McQueen and his group [60,59], developing an improved version for general regions called GENSMAC. An adaptation for generalized Newtonian flow then followed [89], viscous jet buckling was studied using GENSMAC (see [91]) and apparatus was built to experimentally validate the code [92]. At the turn of the millennium saw the development of a solid modeling shell for enhanced visual input and output (see [11]). An axisymmetric version [93] permitted the simulation of G.I. Taylor’s experiments on jets impinging on a quiescent fluid, a full three-dimensional version was produced [11,90,94] and high order monotone upwinding techniques were introduced to permit the simulation of a hydraulic jump (see [23]). In particular the code was shown capable of simulating effects such as the double roller effect (two vortices at the transition layer) predicted by the experiments of Craik et al. [13]. More recently, the GENSMAC code has been extended to cope with generalized Newtonian flows in both two and three dimensions [89,96]. An effective way of dealing with surface tension or interfacial tension has been developed for a 3D version of the code [81]: essentially this involves an accurate computation of both the free surface normal and curvature at a (reasonably large) numbers of points. It has been demonstrated that in practice the code can resolve moving contact angles [57]. However, where the code may prove to be at its most effective is in free surface rheological flows: an Oldroyd-B fluid model has been simulated [95] and several classical problems have been solved.

910

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

3. A review of interfacial tracking techniques The original MAC method was solved on an Eulerian grid and the free surface was defined by whether a cell contained fluid or did not; it was said to contain fluid if the cell contained one or more marker particles and it was defined to be empty if it had no marker particles. Marker and cell methods are now much more sophisticated and provide accurate resolution of interfaces (and free surface) by accurately determining the surface normal and curvature (see [81]). It is crucial, especially when interfacial tension is present, that any algorithm should be capable of not only accurately obtaining and following an interface, but also able to cope with merging and separation. It will be convenient to divide the growing number of methods into four, somewhat overlapping, categories: front tracking, volume of fluid, level set methods and hybrid methods. Shock-capturing methods are usually associated with compressible fluids; these methods are now extremely sophisticated, explicitly enforcing monotonicity through a non-linear step while simultaneously maintaining high order. The reader is referred to, for example, the books of Hirsch [40] and LeVeque [53]. The implementation and development of those ideas are usually credited to Glimm and his co-workers (see e.g. [27] and more recent papers [28–31]). Essentially, they represent the moving front by a connected set of points, which form a moving internal boundary. To calculate the evolution inside the fluid, in the vicinity of the interface, an irregular grid is constructed and a special finite difference stencil is used on these irregular grids. Authors who have used this approach for different flow regimes include Chenn et al. [12], Darip et al. [17], Moretti [63], Peskin [70] (see also Fauci and Peskin [22] and Fogelson and Peskin [25]). Tryggvason and his co-workers [20,34] have used front tracking methods for capturing the deformation of bubbles. More recently, Galaktionov et al. [26] introduced an adaptive technique, based on both surface stretching and surface curvature analysis, for tracking strongly deforming fluid volumes. Shin and Juric [78] employed the constraint of a level contour for front tracking without connectivity. With echoes of Los Alamos, Jaeger and Carin [46] presented a front tracking ALE method in order to compute discontinuous solutions on unstructured finite element meshes. Using an arbitrary Lagrangian–Eulerian formulation, the mesh is continuously adapted by moving the nearest nodes to the interface providing a sharp solution there with no smearing. Rather than employing virtual marker particles, the volume of fluid (VOF) approach (see e.g. [44,52,65,98]) involves using the volume fraction F ðw; tÞ to represent the free surface or interface. Typically, it represents a computational cell wij , e.g. wij ¼ fxi ; xi 6 x 6 xiþ1 g. If F ðw; tÞ ¼ 0 it is void of fluid. If 0 6 F ðw; tÞ 6 1 then the w contains some fluid. An advantage of representing the free surface as a volume fraction is the fact that we can write accurate algorithms for advecting the volume fraction so that mass is conserved, while still maintaining a reasonably sharp

representation of the free surface. On the other hand, a disadvantage of the volume of fluid approach is the fact that it is difficult to accurately compute local curvature from volume fractions. Some authors (e.g. [82,84]) have attempted to overcome this difficulty, but it still remains a problem. Li et al. [54] have used the VOF approach to consider a bubble subjected to shearing between two plates. A modification of the VOF scheme was introduced by Harvie and Fletcher [45], while Zhao et al. [107] developed a VOF approach based on high order characteristics and designed to be employed on unstructured grids. More recently, VOF methods have been developed by Lorstad and Fuchs [55,56] and Kim and Lee [48,49]. Quadtree adaptive VOF methods are also finding favour (see [33,103]). Adaptive grid local refinement has also been employed to increase the accuracy of the VOF method (see [87]). Another approach that has been popular in recent years is the level set method. In this method, the interface is simply captured by a contour of a particular scalar function. This was originally introduced by Osher and Sethian [69], and Sethian has produced two useful books on the subject [76,77]. To date, this interface technique has been used to model interfacial phenomena in fluid dynamics, grid generation and material engineering, to name a few. In the level set method a smooth function, called the level set function, is used to represent the free surface. Liquid regions are regions which have uðx; tÞ > 0. The interface (or free surface) is the set of points such that uðx; tÞ ¼ 0. One of the advantages of the level set method is its simplicity, especially when computing the curvature of a moving free surface. By representing the free surface in this way, the unit normal vector n and the mean curvature j are simply n¼

ru jruj

ð1Þ

and j¼r

ru ; jruj

ð2Þ

respectively. The level set method has the principal disadvantage that the discretization of ut þ u  ru ¼ 0

ð3Þ

(where u ¼ ½uðx; y; z; tÞ; vðx; y; z; tÞ; wðx; y; z; tÞ is the underlying velocity field) is likely to incur more numerical errors than a front tracking or a volume of fluid method, particularly if the free surface experiences rapid changes in curvature. A common problem is that mass is not conserved. Recent papers on level sets applied to free surface and interfacial problems include Sussman and Fatemi [83], Quecedo and Pastor [73] (see also [71]), Yue et al. [106] and Kohno and Tanahashi [50]. Level set methods cope well with merging and separation of interfaces. They provide a direct means of obtaining the normal and curvature, but they suffer from inaccuracy and smearing. On the other hand, VOF, particle methods and front tracking, while able

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

to give, generally speaking, reasonably accurate representation of the interface, cannot easily obtain the normal and the curvature locally. It is therefore perhaps not surprising that a number of hybrid approaches have appeared recently. A level set approach using Lagrangian marker particles has been suggested by Enright et al. [19]. A VOF approach using ‘‘baby” cells have been put forward by Kim and Lee [48]. A new VOF scheme that preserves mass exactly has been developed by Aulisa et al. [2] (see also [33,47,55]). Nourgaliev et al. [66] introduce a characteristics-based matching method which capitalizes on recent developments on level sets. A novel scheme has been developed by coupling level sets with adaptive mesh refinement (see [50]), while [3] combines a front tracking approach with marker particles which are area preserving. The boundary integral method is another approach that may be used to resolve the free surface and can be effective when inertial forces are negligible. In this method, the flow equations are mapped from the immiscible fluid domains to the sharp interfaces separating them, thus reducing the dimensionality of the problem (the computational mesh discretizes only the interface). In particular, Baker and Moore [4] and Tsai and Miksis [101] have used this approach to solve axisymmetric bubbles rising in a liquid. Impressive results were also obtained for the problem of bursting bubbles by Boulton and Blake (see [6]). 4. The emphasis of the review

911

ru¼0

ð4Þ

and conservation of momentum     Du ou ou q ¼q þ u  ru ¼ q þ r:ðuuÞ Dt ot ot ¼ r  r þ qg ¼ rp þ r  s þ qg;

ð5Þ

where u is the fluid velocity field and r is the symmetric stress tensor; p and q are, respectively, the pressure and density of the fluid, and qg is the gravitational force. The tensor s is known as the extra-stress tensor. If multiphase flow is being considered then both density and viscosity are deemed to be variable over the whole domain, but constant in each specific fluid. Thus each fluid is taken to be incompressible and the continuity Eq. (4) is valid over the whole domain. The pressure is discontinuous along the interfaces, with a discontinuity proportional to the local surface tension r and curvature j. Hence it behaves like a Heaviside function and consequently the term rjrH

ð6Þ

needs to be added to the right-hand side of Eq. (5). Thus the interface(s) can be identified where the gradient of this Heaviside function is non-zero. A reasonably general constitutive relationship is that for an Oldroyd-B fluid   O O s þ k1 s ¼ 2l D þ k2 D ; ð7Þ O

Having provided a historical backdrop to the MAC method, the purpose of this paper is to describe the essence of the modern marker and cell approach in sufficient detail to convey the flavor of the method, while not overburdening the reader with all the minor details associated with the methodology. For example, in three dimensions, to approximate the free surface in a reasonably accurate manner, it is necessary to consider 24 cases; in this review we only discuss this in generality, referring the reader to the appropriate paper. The same is also true of the many approximations to curved boundaries. Equally, there are many bounded monotone upwinding schemes available; we focus on the variable-order non-oscillatory scheme (VONOS) [102] and refer the reader to suitable references. This paper neither contains code verification nor model validation, but refers the reader to previous papers. However, the review does include a number of applications, by no means exhaustive, to illustrate the power and potential of the method. The paper concludes with a speculative discussion on the future of this and related approaches to free and interfacial flows and proposes six test problems that a good free surface flow code should be capable of solving. 5. Governing equations The basic equations governing the flow of an incompressible fluid are the conservation of mass

where the upper convected derivative A is defined by O

A¼

DA T  ðruÞ  A  A  ðruÞ: Dt

The rate of deformation tensor is  1 D ¼ ðruÞ þ ðruÞT 2

ð8Þ

ð9Þ

and k1 and k2 are time constants (relaxation and retardation respectively) and l is the apparent viscosity (note if l is constant it is the dynamic viscosity). It will be found useful to decompose the extra-stress tensor into Newtonian and non-Newtonian parts   k2 s ¼ 2l D þ S; ð10Þ k1 where S represents the non-Newtonian contribution of the extra-stress tensor. If in Eq. (7) k2 ¼ 0 we obtain the Maxwell model (pure elasticity). If we choose k1 ¼ k2 ¼ 1 we obtain the Navier– Stokes equations of Newtonian fluid flow. The generalized Newtonian flow equations arise if, in Eq. (7), k1 ¼ k2 ¼ 0 but l is now a function of the local shear rate q. Shear thinning fluids occur frequently in industrial processes and this model is often sufficient to provide realistic solutions. If qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q ¼ 2trðD2 Þ with D given by (9), then we may employ the Cross model [14]

912

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

l  l1 1 ¼ m; 1 þ ðKqÞ l  l0

ð11Þ

where m, l1 , l0 and K are given positive constants related to the particular fluid under study. Note also that other alternatives to the Cross model include Carreau–Yasuda and power-law models [5]. Note also that in the rheological literature it is more common to use c_ (rather than q) for the local shear rate. Non-dimensionalising with respect to appropriate length (L) and velocity (U) scales, with q0 , l0 and r0 appropriate reference values for the density, viscosity and surface (interfacial) tension respectively, leaves the continuity equation unchanged, but transforms (5) into Du 1 ¼ rp þ r  s þ g; Dt Fr

ð12Þ



ou 1 1 þ r:ðuuÞ ¼  rp þ r  l ru þ ðruÞT ot q qRe 1 j þ 2 g  rH xe Fr

ð20Þ

with now Re ¼ qUL=l. 5.1. Boundary conditions The boundary conditions at the mesh boundary can be of several types, namely, no-slip, free-slip, prescribed inflow, prescribed outflow and continuative outflow. Let un , ul and um denote the normal and tangential velocities to the boundary, respectively. Then, for a no-slip rigid boundary we have un ¼ 0;

ul ¼ 0;

um ¼ 0

pffiffiffiffiffiffiffiffiffiffiffi where Fr is the Froude number given by Fr ¼ U = 2Ljgj. The surface tension (6) becomes

and for a prescribed inflow

1 rj rH ; xe q

where the subscript n, l and m denote the normal and the two tangential directions to the boundary, respectively. Of course, if the solid boundary had both a rotational and/or a translational velocity, then this could equally be described by no-slip conditions. For the Poisson equation (see Section 7) we require, for the potential function w, homogeneous Dirichlet boundary conditions at the free surface and homogeneous Neumann boundary conditions at the solid boundaries.

ð13Þ

where xe is the Weber number which is xe ¼

q0 LU 2 ; r0

while (7) becomes   5 5 2l s þ We s ¼ D þ k2 D ; qRe

ð14Þ

ð15Þ

un ¼ U inf ;

ul ¼ 0;

um ¼ 0;

5.2. Interfacial and free surface conditions

where We ¼ UL1 k1 is the Weissenberg number and Re ¼ q0 ULl1 0 is the Reynolds number. Note that in order to avoid possible confusion between the Weber and Weissenberg numbers a different font has been used. Eq. (9) remains unchanged, but (10) becomes   2l 1 k2 s¼ D þ S: ð16Þ q Re k1

As we shall see, the proposed method is capable of dealing with a low density fluid in two different ways: as one of the phases of the flow or as an inert phase, thus introducing a moving free boundary. At the free surface, the boundary conditions for pressure and velocity, assuming zero viscous stress in the empty phase, are given by

Of course when there is only one phase, the density and viscosity are constant and the surface tension is either constant or zero, we may choose q0 ¼ l0 ¼ r0 ¼ 1 and (13) becomes

where n, l and m are the local normal and tangential vectors to the free surface/interface. The viscous stress–tensor in the Newtonian case is given by

1 jrH xe with now xe ¼ qLU 2 =r and (16) becomes   2 k2 s¼ D þ S: Re k1 The constitutive Eq. (15) now becomes   5 k2 S þ We S ¼ 2 1  D: k1

ð17Þ

ðr  nÞ  n ¼ P cap ;

r ¼ pI þ

2l D Re

ðr  nÞ  l ¼ 0;

ðr  nÞ  m ¼ 0;

ð21Þ

ð22Þ

and P cap ¼ rjWe1 is the capillary pressure, originating from the effects of surface tension r. 6. Discretization

ð18Þ

ð19Þ

If further S  0, then substitution of (18) into (5) results in the familiar Navier–Stokes equations

The MAC method employs finite differences on a staggered grid. For non-Newtonian fluids this can become extremely complex and so we shall focus on the discretization of the Navier–Stokes Eq. (20). Fig. 1 displays an example of a staggered grid cell and the position of the variables in this cell: the velocities are calculated at the faces, and the other variables (pressure, density, viscosity

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

913

cell, and we compute an approximation by averaging the known neighbouring values. Thus 2ðpiþ1;j;k  pi;j;k Þ 1 op ¼ : q ox Dxðqiþ1;j;k þ qi;j;k Þ

ð28Þ

To compute the viscosity, we use the harmonic mean ! 1 1 : ð29Þ liþ12;j;k ¼ 2 þ liþ1;j;k li;j;k This is more accurate than a simple average [100]. Thus, the viscous term V 1x using (29) becomes Fig. 1. Typical cell showing the location of the dependent variables.

 V 1x iþ1;j;k ¼ 2

and, in the case of non-Newtonian flow, the non-Newtonian part of the extra-stress tensor – in this figure these variables are denoted by /) are computed at the cell centre. Consider the x component of the non-dimensional momentum equation, expressed in Cartesian coordinates as ou oðuuÞ oðuvÞ oðuwÞ þ þ þ ot ox oy oz  2  1 op l o u o2 u o 2 u þ ¼ þ þ q ox qRe ox2 oy 2 oz2      1 ol ou ol ov ou ol ow ou 2 þ þ þ þ þ qRe ox ox oy ox oy oz ox oz 1 rj oH : þ 2 gx  qxe ox Fr

2

ð24Þ

where oðuuÞ oðuvÞ oðuwÞ þ þ ; ox oy oz   l o2 u o2 u o2 u V 1x ¼ þ þ qRe ox2 oy 2 oz2

1

ð30Þ

 V 2x iþ1;j;k ¼

This equation can be written in the compact form

Cx ¼

4



Re liþ1;j;k þ li;j;k qiþ1;j;k þ qi;j;k  ui12;j;k  2uiþ12;j;k þ uiþ32;j;k  Dx2 uiþ1;j1;k  2uiþ12;j;k þ uiþ12;jþ1;k þ 2 Dy 2  uiþ1;j;k1  2uiþ12;j;k þ uiþ12;j;kþ1 þ 2 : Dz2 1

The value of the viscosity on the edges of the cell is computed in a similar manner to (29). Thus V 2x becomes

ð23Þ

ou 1 op 1 rj oH þ Cx ¼  þ V 1x þ V 2x þ 2 gx  ; ot q ox qxe ox Fr



ð25Þ ð26Þ

and      1 ol ou ol ov ou ol ow ou 2 þ þ þ V ¼ þ qRe ox ox oy ox oy oz ox oz 2 x

ð27Þ are the convective (C x ) and viscous (V 1x , V 2x ) terms of (23). Eq. (23) is discretized by finite differences at the point ði þ 12 ; j; kÞ of the cell as follows. The time derivative is discretized by the first-order forward difference. The other terms of (24) are discretized at the time level n, i.e. the discretization is explicit. Thus, in the following equations, the index n is omitted to simplify the notation. The pressure gradient is discretized by central differences. Notice that it is necessary to evaluate q at the face ði þ 12 ; j; kÞ of the

1 Reðqiþ1;j;k þ qi;j;k Þ ( ðliþ1;j;k  li;j;k Þðuiþ32;j;k  ui12;j;k Þ  Dx2 2 !1 4 4 1 1 1 1 þ þ þ þ Dy li;j;k liþ1;j;k li;jþ1;k liþ1;jþ1;k !1 3 1 1 1 1 5  þ þ þ li;j;k liþ1;j;k li;j1;k liþ1;j1;k  viþ1;jþ12;k þ viþ1;j12;k  vi;jþ12;k  vi;j12;k  2Dx  uiþ12;jþ1;k  uiþ12;j1;k þ 2Dy 2 !1 4 4 1 1 1 1 þ þ þ þ Dz li;j;k liþ1;j;k li;j;kþ1 liþ1;j;kþ1 !1 3 1 1 1 1 5  þ þ þ li;j;k liþ1;j;k li;j;k1 liþ1;j;k1  wiþ1;j;kþ12 þ wiþ1;j;k12  wi;j;kþ12  wi;j;k12  2Dx ) uiþ12;j;kþ1  uiþ12;j;k1 þ : ð31Þ 2Dz

The discretizations of the momentum equations in the y and z directions are obtained in a similar fashion. The non-linear terms in the momentum equation are discretized

914

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

using the high order bounded upwind VONOS scheme, omitted here for brevity. Full details of the implementation may be found in Ferreira et al. [23]. This scheme is accurate and stable for simulations where the Reynolds number is high. 6.1. Markers and cells The essence of the MAC method is the virtual marker particles and the cells defined on an Eulerian grid. Marker particles, often as many as 16 per cell, are moved from their coordinates at time t to their new coordinates at time t þ Dt according to the newly computed velocity u at the cell centre. The flow particles are evaluated on a three-dimensional uniform Eulerian grid, in which every cell, at each time step, is classified according to its position relative to the fluids and the rigid, or indeed moving, boundaries or/and interfaces. Cells with more than half their volume within the computational domain are classified as BOUNDARY (B) cells; the same criterion is used for INFLOW (I) cells. Any cell completely inside the fluid is classified as a FULL (F) cell, those completely outside the fluid (but still within the computational domain) are EMPTY (E) cells. Those containing marker particles are SURFACE (S) cells if they have at least one E neighbour. These criteria are applied to each fluid in the simulation, and the INTERFACE cells are identified as the cells that are S cells or F cells for more than one fluid at the same time. Fig. 2 shows an example of cell flagging in a two-dimensional case (simply for clarity), with a classification for each fluid. The cell classification is updated at each time step using information provide by the (virtual) Lagrangian mesh constituted by the marker particles. Since these particles are restricted by a CFL type condition, the interfaces move at most one cell at each time step. Hence, it is only necessary to update cell classification in the neighbourhood of the surface cells. The procedure can be described in three steps. First, the cells that contain particles are identified. If any of these identified cells are, at the next time step, found to be empty (E) then they are marked as surface (S) cell. Secondly, cells previously marked as S that no longer contain particles are reflagged as either E or F

according to the status of their neighbours. Thirdly, the full cells that contain particles and have any E neighbouring cells are marked as S. (A neighbouring cell is defined as any cell that is contiguous with the face of another cell; cells that touch at edges or vertices are not considered to be neighbours.) 6.2. Initial approximation of the free surface In order to apply the free surface boundary conditions (21) in each S cell, approximations for the surface normals are required. These are obtained according to the classification of the neighbouring cells (see e.g. [81]). For example: nc ¼ ð1; 0; 0Þ, nc ¼ ð0; 1; 0Þ or nc ¼ ð0; 0; 1Þ if only one neighbour is an E cell; nc ¼ ð p1ffiffi2 ;  p1ffiffi2 ; 0Þ, nc ¼ ð0;  p1ffiffi2 ;  p1ffiffi2Þ or nc ¼ ð p1ffiffi2 ; 0;  p1ffiffi2Þ if there are two neighbouring cells; and nc ¼ ð p1ffiffi3 ; 0; 0Þ, nc ¼ ð0;  p1ffiffi3 ; 0Þ or nc ¼ ð p1ffiffi3 ;  p1ffiffi3 ;  p1ffiffi3Þ if there are three neighbouring cells. Note that the notation nc means the normal at the centre of the S cell under consideration. When surface or interfacial tension is absent this has proven to be an effective approach. However, when it is present it is not only necessary to obtain more accurate estimates of the normals: it is necessary to obtain a good approximate of the surface curvature at the centre of each surface cell. To do this, it is crucial that sub-cell surface/interfacial tension effects are taken into account. In Section 9 we shall show how this is achieved and how it can be used in the computation of the capillary pressure. 6.3. Interfacial tension effects In order to implement the interfacial tension effects in the momentum equations, it is necessary to determine the resulting faces that need to be applied in the discretized domain. If the mesh is chosen to be sufficiently fine, then the interface can be approximated by the cells containing the front. The interfacial tension forces in those cells are given by a discretized form of rj f¼ rH ; ð32Þ qxe

Fig. 2. Multi-fluid cell flagging example: (a) classification for fluids 1 and 2; (b) classification only for fluid 1; (c) classification only for fluid 2.

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

where H is a Heaviside function that is 1 inside the particular fluid under consideration and 0 outside. It is convenient to define one Heaviside function H k for each fluid k in the discretized domain. The Heaviside function is constructed at each time step according to the cell classification, and is given by 8 > < 1 if ðx; y; zÞ is inside a full cell or H k ðx; y; zÞ ¼ an interface cell of the fluid k; > : 0 otherwise;

915

oH 2 H 2 ðiDx; jDy; kDzÞ  H 2 ðði  1ÞDx; jDy; kDzÞ j 1 ¼ Dx ox ði2;j;kÞ 1 ¼ Dx ð35Þ and then the resulting force to be applied to face ði  12 ; j; kÞ is given by !   r oH 1  oH 2  j1 þ j2 fi12;j;k ¼  2qi12;j;k We ox ði1;j;kÞ ox ði1;j;kÞ 2

such that the force to be applied in each cell is given by f¼

n r 1X jk rH k : qxe n k¼1

ð33Þ

Each fluid supplies a different curvature and Heaviside function. For each fluid k the force is applied at each face of the interfacial cell. For example, to calculate the force contribution in a single face ði  12 ; j; kÞ, as shown in Fig. 3, we have two Heaviside functions to be discretized, H 1 and H 2 . oH 1 H 1 ðiDx; jDy; kDzÞ  H 1 ðði  1ÞDx; jDy; kDzÞ j 1 ¼ Dx ox ði2;j;kÞ ¼ 0; ð34Þ

rj2 ¼ : 2Dxqi12;j;k xe

ð36Þ

The other cases follow a similar procedure. It is the gradient of the Heaviside functions that provides the correct direction of the applied forces in any particular configuration. 6.4. Surface merging and splitting We describe a simple technique which works well in the absence of surface tension. The ideas will be presented pictorially. 6.4.1. Merging Consider the example of two free surfaces traveling towards each other as depicted in Fig. 4a. Without loss

Fig. 3. Interface cell between fluids 1 and 2, with full cells in the opposite faces.

a

2

b

Fig. 4. Steps required to perform surface merging.

c

916

a

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

b

c

Fig. 5. Steps required to perform surface splitting.

of generality assume that the left free surface is traveling faster than the one on the right. Initially at time tn they are still well separated. However, at time tnþ1 two full cells associated with the left free surface have become surface cells, resulting in both free surfaces having surface cells in common. When this occurs the two free surfaces are deemed to have merged, the common surface cells are changed to full cells (see Fig. 4b) and then the free surface is redrawn (see Fig. 4c). Note at this point the new free surface(s) is non-smooth. However, upon the next calculational time step, the marker particles will move creating a smooth free surface again. 6.4.2. Demerging We shall illustrate demerging or splitting with the neck of a pendant drop. As the neck becomes thinner it will reach a point where there may be only two full cells within (the narrowest part of) the neck (see Fig. 5a). After the next calculational cycle (see Fig. 5b) (or indeed several calculational cycles) this may reduce to two contiguous surface cells. When this occurs the fluid is deemed to have demerged, that is the neck has split. The single free surface of the pendant drop divides into two free surfaces as shown in Fig. 5c. Again momentarily the free surfaces are nonsmooth, but this is resolved after the next calculational cycle. When surface tension is present and dominant, this approach may prove inadequate. How this is best implemented in three-dimensional flows is still an open question. However, progress has been made with other front tracking methods and the reader is referred to [100].

Step 1: Let ~pðx; tÞ be a pressure field satisfying the correct pressure condition on the free surface. This pressure field is computed according to the required boundary stress conditions. Step 2: The intermediate velocity field ~uðx; tÞ is computed by the explicit1 discretized form of the momentum equations o~u 1 1 þ r:ðuuÞ ¼  r~p þ r  2lD ot q qRe 1 rj þ 2g rH qxe Fr

with ~uðx; t0 Þ ¼ uðx; t0 Þ using the correct boundary conditions for uðx; t0 Þ. It can be shown [94] that ~uðx; tÞ possesses the correct vorticity at time t. However, ~uðx; tÞ does not satisfy (4) in general. Let 1 uðx; tÞ ¼ ~uðx; tÞ  rwðx; tÞ q with   1 rwðx; tÞ ¼ r  ~uðx; tÞ: r q

ð38Þ

ð39Þ

Thus uðx; tÞ satisfies (4) and the vorticity remains unchanged. Therefore, uðx; tÞ is identified as the updated velocity field at time t. Step 3: Solve the elliptic Eq. (39). Step 4: Compute the velocity field (38). Step 5: Compute the pressure using

7. Numerical procedure pðx; tÞ ¼ ~pðx; t0 Þ þ The numerical procedure is based on the GENSMAC3D method [94], but involves the essential ideas of the original SMAC algorithm [36]. It is supposed that the velocity field uðx; t0 Þ is known at a given time t0 , and that the boundary conditions for the velocity and pressure are given. The updated field uðx; tÞ at t ¼ t0 þ Dt is calculated using the following algorithm:

ð37Þ

wðx; tÞ : Dt

ð40Þ

Step 6: Update the position of the marker particles by solving x_ ¼ uðx; tÞ by using Euler’s method, or similar explicit scheme. 1 Of course this is not necessary and an implicit formulation can lead to considerable computational savings (see [67]; see also [18,72,32]).

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

The last step in the calculation involves moving the marker particles to their new positions. The coordinates of these virtual particles are stored and updated at each time step by solving the differential equation dx=dt ¼ uðx; tÞ. This provides the new coordinates of every particle, allowing us to determine whether or not it has moved to a new computational cell or it has left the containment region through an outlet. By using the front tracking methodology [81], only marker particles on the free surface and the interfaces need to be considered. The elliptic Eq. (39) is solved using the conjugate gradient method with the previous velocity field employed as the initial guess. However, as the density variation across the interface increases, more iterations are required for convergence of the method. In these cases, it is advisable to use a diagonal preconditioner to speed up the convergence of the conjugate gradient method. The numerical procedure we have just described implicitly assumes that we are dealing with Newtonian fluids. The algorithm extends fairly naturally to non-Newtonian fluids (and generalized Newtonian fluids, see [96]). Two further steps are required between step 4 an step 6, namely: Step 5a: Update the components of the non-Newtonian extra-stress tensor on the rigid boundaries (for more detail, see [96]). Step 5b: Compute the non-Newtonian extra-stress tensor S from Eq. (19)   O k2 S þ We S ¼ 2 1  D: ð41Þ k1 Eq. (41) is solved by finite differences in an analogous way to the Navier–Stokes equation. Details will not be given here, but it can be found in [92] for the two-dimensional case and in Tome´ et al. [97] for the three-dimensional case.

7.1. Time-stepping procedure A time-stepping procedure for computing the appropriate time-step size for every cycle is employed. It is based on the stability conditions (written in non-dimensional form) Dt < Dt <

Dx ; kuk Dx2 Dy 2

ð42Þ Dx2 Dy 2 Dz2 Re ; þ Dx2 Dz2 þ Dy 2 Dz2 2

ð43Þ

where the first inequality is understood componentwise. The restriction (42) requires that no particles should cross more than one cell boundary in a given time interval; this is an accuracy requirement. The second restriction (43) comes from the explicit discretization of the Navier–Stokes equations and is essentially a local von Neumann stability requirement. Since low Reynolds number flows (0 6 Re 6 10) are the primary concern, it is anticipated that (43) is generally the more restrictive condition.

917

Let Dtn be the time-step employed in the previous calculational cycle and define 1 Dx2 Dy 2 Dz2 Re ; 2 Dx2 Dy 2 þ Dx2 Dz2 þ Dy 2 Dz2 2 1 Dx Dtu ¼ k1   ; 2 je u max j 1 Dy Dtv ¼ k2   ; 2 jev max j 1 Dz Dtw ¼ k3   ; e max j 2 jw Dtvisc ¼

ð44Þ ð45Þ ð46Þ ð47Þ

where j~umax j, j~vmax j and j~ wmax j are maximum of the tilde velocities computed through (37) and 0 < ki 6 1; i ¼ 1; 2; 3. The extra factor of 0.5 in (45)–(47) has been introduced as a conservative measure to allow for the fact that only local stability analyses have been performed. The time-step employed in the calculation is then chosen to be Dtnþ1 ¼ k  minfDtvisc ; Dtu ; Dtv ; Dtw g;

ð48Þ

where 0 < k 6 1. The factor k in (48) is necessary since the values of junþ1 max j, nþ1 nþ1 jvmax j and jwmax j at the beginning of the calculational cycle are not known. It is computationally more efficient to use the tilde velocities with the factor k acting as a compensatnþ1 nþ1 ing measure than using umax , vnþ1 max , and wmax . However, if n Dtu , Dtv or Dtw is less than Dt then the tilde velocities are recalculated and the time-step is revised. This time-stepping procedure has proved to be an efficient automatic means of adjusting the size of the time-step. 7.2. Finite difference approximation of the convective terms For higher Reynolds numbers (Re  1000) the convective terms must be approximated by monotone, bounded, high accuracy upwinding. The upwinding scheme we have selected to use is that due to Varonos and Bergeles [102] (known as VONOS). There are of course many others, but this proved to be one of the most effective when solving the hydraulic jump problem (see Section 12.3) – the reader is referred to the comparison made by Ferreira et al. [23]. We require an approximation to the convective term, r  ðuuÞ, that is for the x component a bounded monotone high accuracy finite difference approximation is required for    oðuuÞ oðuvÞ oðuwÞ þ þ : ox  1 oy  1 oz  1 iþ2;j;k

iþ2;j;k

iþ2;j;k

Since the VONOS scheme is complex and lengthy, it is omitted here; the reader may wish to consult the original paper [102] or Ferreira et al. [23] for details. 8. Undulation removal The computation of the surface and interfacial tension is performed at two levels: first at the sub-grid level, where

918

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

small undulations in the fronts (free surfaces and interfaces) are eliminated; and second at the cell level where the curvature at each surface or interface cell is approximated. This approximation will be used in the implementation of the pressure boundary condition at the free surface and in the calculation of the interfacial forces discussed in Section 8.3. In many applications, in particular when the Reynolds number is high (larger than 50), small undulations may appear at the free surface and at the interfaces. This is due to variations in the velocity field from cell to cell and may be amplified in regions where the surface area is shrinking. They are physically removed by a combination of surface and interfacial tension and viscous effects. A numerical implementation that acts at cell level cannot take into account these sub-cell undulations and correctly suppress then. There are several approaches that can be used to suppress these unphysical undulations such as applying a Gaussian filter. However, in fluid flow simulations it is vital that, whatever technique is used, it does not alter the mass (or, equivalently, the volume in the case of an incompressible fluid). The filter we describe, TSUR-3D (TSUR: Trapezoidal Sub-grid Undulation Removal) has been designed to conserve volume (mass). It is applied to an unstructured grid, represented by a ‘‘half-edge” data structure [11]. Fig. 6 shows a typical vertex v and its corresponding star formed by the set of vertices xi ; i ¼ 0; 1; . . . ; n, which are connected to it with typical edges ðv; xi Þ and ðv; xiþ1 Þ and typical faces ðxi ; v; xiþ1 Þ, i ¼ 0; 1; . . . ; n, where xnþ1 is to be interpreted as x0 . Let p be an arbitrary point. Then the volume of the star with vertices v and p is the volume bounded by the faces ðxi ; v; xiþ1 Þ, i ¼ 0; 1; . . . ; n, respectively.

8.1. Normal vector calculation Both vertex balance and undulation removal procedures require a normal direction which needs to be determined from the neighbouring vertices xi , i ¼ 0; 1; . . . ; n. Let v be a vertex of the interface and let SðvÞ be the star of v. Let p be the average of the vertices in SðvÞ given by n 1X p¼ xi ; ð49Þ n i¼0 where x0 ; x1 ; . . . ; xn 2 SðvÞ. Thus we can compute the normal vector by the average of the normal vector at the faces ðxi ; p; xiþ1 Þ, i ¼ 0; 1; . . . ; n, given by Pn1 ðxi  pÞ ^ ðxiþ1  pÞ þ ðxn  pÞ ^ ðx0  pÞ n ¼ Pi¼0 ; n1 k i¼0 ðxi  pÞ ^ ðxiþ1  pÞ þ ðxn  pÞ ^ ðx0  pÞk ð50Þ where ^ denotes the vector (cross) product. Note this is a weighted average: normals of those faces with greater area will be given greater weight. This equation provides a robust normal direction for this method. 8.2. Vertex balance procedure The vertex balance procedure is applied at each vertex of the fluid in order to move that vertex to some new position on a line passing through p given by (49) with the direction n given by (50), such that the local volume is preserved. Let p be a plane that contains p with normal vector n. If p is inside the projection of SðvÞ onto p, we compute the new position of v by vnew ¼ p þ hn;

ð51Þ

where h is determined such that the local volume is preserved. This is guaranteed if V ¼ hV 1

or

h¼

V ; V1

ð52Þ

where V is the volume of the polyhedron ðv; x0 ; x1 ; . . . ; xn ; pÞ and V 1 is the volume of the unitary polyhedron ðp þ n; x0 ; x1 ; . . . ; xn ; pÞ. This follows since the volume of the polyhedron ðp þ hn; x0 ; x1 ; . . . ; xn ; pÞ is equal to h times the volume of the unitary polyhedron ðp þ n; x0 ; x1 ; . . . ; xn ; pÞ. The volume of the tetrahedron ðx1 ; x2 ; x3 ; x4 Þ is given by 2 3 ðx2  x1 Þ ðx3  x1 Þ ðx4  x1 Þ 1 6 7 ð53Þ V ¼ det 4 ðy 2  y 1 Þ ðy 3  y 1 Þ ðy 4  y 1 Þ 5; 6 ðz2  z1 Þ ðz3  z1 Þ ðz4  z1 Þ

Fig. 6. A typical vertex and its star projected into the plane p.

where xi ¼ ðxi ; y i ; zi Þ; i ¼ 1; . . . ; 4. Since the volume of any polyhedron may be considered to be the sum of the volume of several tetrahedra, this can be straightforwardly calculated. As we can see in Fig. 7 this procedure does not remove any undulations, but just moves the vertex to a new position that is better centred in its star, improving the quality

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

919

Thus, we have two polyhedra ðv1 ; x10 ; x11 ; . . . ; x1n ; p1 Þ with volume V 1 and ðv2 ; x20 ; x21 ; . . . ; x2n ; p2 Þ with volume V 2 , where ðx10 ; x11 ; . . . ; x1n Þ 2 Sðv1 Þ and ðx20 ; x21 ; . . . ; x2n Þ 2 Sðv2 Þ (see Fig. 8). We have to determine the new height h in the direction of n such that the local volume is preserved. We can exploit linearity between volumes of polyhedra and their heights by writing V 1 þ V 2 ¼ V 1 ðh=h1 Þ þ V 2 ðh=h2 Þ

ð56Þ

or Fig. 7. A vertex with its star, before and after the vertex balance procedure has been applied.

h¼

V1þV2 : V 1 =h1 þ V 2 =h2

ð57Þ

The new positions of the vertices v1 and v2 are given by of the surface mesh while preserving local volume. As we shall see this improves the robustness of the undulation removal procedure presented in the next section. 8.3. Undulation removal procedure The vertex balance procedure is designed to produce a more homogeneous unstructured grid that represents the free surface or interface. However, it cannot remove the sub-grid undulations; a smoothing procedure is required for the effective removal of these undulations. Let e be an edge of a surface and v1 and v2 its vertices as shown in Fig. 8. The smoothing procedure changes the positions of v1 and v2 simultaneously in a such way that the local volume is preserved. Let n1 and n2 be the normal vectors computed by (50) at the vertices v1 and v2 , respectively. We consider a normal vector at the edge e given by the average of n1 and n2 , namely: n¼

n1 þ n2 : kn1 þ n2 k2

h1 ¼ hv1  m; ni and

h2 ¼ hv2  m; ni;

ð54Þ

where h; i denotes the inner (or dot) product. Next, we compute the points p1 and p2 by and

p2 ¼ v2  h2 n:

ð55Þ

ð58Þ

This undulation removal procedure is applied periodically to all the edges of the surface mesh with a period that has to be chosen such that the surface is sufficiently smooth, and the large scale undulations are not greatly affected. The frequency of application is dependent on the particular problem, the grid resolution and the time step. Of course,this undulation removal procedure and the curvature calculation of the next section need not be applied to a large class of problems – non-Newtonian low Reynolds number flows, especially when surface or interface tension can be neglected. 9. Curvature calculation The curvature is approximated by the surface that best fits the points (the positions of the marker particles) in that cell and its neighbours. Numerical experimentation (see [80,81]) suggests approximation by the paraboloid pðx; yÞ ¼ ax2 þ bxy þ cy 2 þ dx þ ey þ f :

Let m ¼ 12 ðv1 þ v2 Þ be the average of the vertices of the edge e. We determine the heights h1 and h2 in the direction of this normal vector by

p1 ¼ v 1  h1 n

v1 ¼ p1 þ hn and v2 ¼ p2 þ hn:

ð59Þ

In this approximation, we first need to determine the normal vector at the cell centre. We can summarize the method of curvature approximation by the following steps: 1. Given a surface at an interface cell, consider all points that lie inside a sphere S d with its centre at the cell centre and a given radius d. 2. Fit a plane to these points and compute its normal. 3. Take a new coordinate system normal to the surface using the computed normal from the second step.

Fig. 8. Undulation removal procedure, showing an edge and its star.

920

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

4. Consider a new sphere S  , centred at the cell centre, with a specified radius  and then map the points in this sphere to the new coordinate system. 5. In the new coordinate system compute the coefficients of the surface given by (59). 6. Compute the curvature of the fitted surface using j¼

^ o2 p on2 p 1 þ ðo^ Þ on

2

3=2 þ

^ o2 p og2 p 1 þ ðo^ Þ og

2

ð60Þ

3=2 ;

^ is the fitted surface given by (59) and ðn; g; fÞ are where p the coordinates of the new coordinate system. Note the f-axis is parallel to the normal vector. In this method, the quality of the approximation is directly related to the chosen radii d and . Typically, d is chosen to be of the order of a cell size and , 1.5 times larger. For more detailed discussion of the influence of these constants see [81], where the optimum choice for these numbers is discussed. How the normal vector and the curvature approximations are made will be described in detail in the next subsection. 9.1. Normal vector approximation The region S d associated with the cell S will include a number of points (coordinates of the marker particles). Let that number be m ¼ mðdÞ. Let xi ¼ ðxi ; y i ; zi Þ, i ¼ 1; 2; . . . ; md , be the number of points inside S d and thus the number of points in cell S and its surrounding neighbours. The equation of the plane will be either z ¼ f ðx; yÞ, y ¼ f ðx; zÞ or x ¼ f ðy; zÞ, according to the maximum absolute value of the cell normal vector components. For example, if jhnC ; ð0; 0; 1Þij > jhnC ; ð1; 0; 0Þij

and

jhnC ; ð0; 0; 1Þij

> jhnC ; ð0; 1; 0Þij; where h; i denotes the inner product, the equation of the plane is given by z ¼ f ðx; yÞ ¼ ax þ by þ c and the coefficients a, b and c are then determined as a best least squares fit to the given points. The normal vector to the plane z ¼ ax þ by þ c can then be calculated from ns ¼

ða; b; 1Þ : kða; b; 1Þk2

ð61Þ

The other cases follow similarly. This procedure determines ns , but the sign may be incorrect. The correct sign can be determined by comparing ns with the normal at the cell centre nc , previously obtained (Section 6.2). If hnC ; nS i > 0 the sign is correct; otherwise the sign must be reversed. 9.2. Curvature approximation Now let xi ¼ ðxi ; y i ; zi Þ, i ¼ 1; 2; . . . ; m , be the points lying in S  . With the normal vector ns given by the plane fit-

ting technique of Section 9.1, we choose another two tangent vectors at the surface such that they form an orthonormal basis together with the normal vector. First we choose two of the three canonical vector e1 ¼ ð1; 0; 0Þ, e2 ¼ ð0; 1; 0Þ and e3 ¼ ð0; 0; 1Þ whose projections on the direction of ns vector are smallest. Second, we apply the Gram Schmidt method to obtain an orthonormal basis. Thus, we have a new coordinate system and the points xi can be expressed in terms of this new coordinate system as ni ¼ ðni ; gi ; fi Þ. This, then, allows us to fit the paraboloid ^ðn; gÞ ¼ an2 þ bng þ cg2 þ dn þ eg þ f p using least squares. This gives us the ‘‘best-fit” values for a, b, c, d, e and f from whence we obtain the curvature ! a b j ¼ 2 þ : ð62Þ 3=2 3=2 ð1 þ d 2 Þ ð1 þ e2 Þ 10. A dynamic contact angle When two fluids, in this case water and air, meet at a solid surface they meet at an angle that is generally not normal to that surface. This is called the contact angle and if the fluids are not quiescent then this angle will move. This Section is concerned with how the MAC method can track this moving contact angle. The contact angle will have the effect of influencing the curvature in the surface cells adjacent to the boundary cells. Thus these surface cells will not be employed directly in the computation of the curvature. Instead the curvature in this case is estimated using the free surface normal n1 , computed at the point x1 of the surface in the adjacent cell opposite to the wall; the coordinates of x1 were previously obtained using the methods of Section 6.1, using the normal to the free surface at the wall (which is prescribed in the case of a constant angle) and the normal distance b of x1 from the wall. The details of the method may become somewhat clearer by referring to the situation depicted in Fig. 9, which displays, for clarity, the equivalent twodimensional problem. Again for simplicity of exposition, we shall, in the following, restrict ourselves to two dimensions. Generalization to three dimensions should be clear. Let S 1 be a surface cell adjacent to a boundary cell B for which we need to compute the capillary pressure. Let n1 be the free surface (interface) unit normal, computed at the centre x1 of the adjacent surface cell S 2 . Let n3 be the wall unit normal, and let n2 be the free surface (interface) unit normal at the contact point which is obtained by adding the (prescribed) contact angle / to the wall angle. Let sw be the unit vector tangent to the wall. Finally, let ss be the unit vector tangent to the free surface at the point of contact x2 at the current time t. The point of contact at a ^2 . previous time ^t ¼ t  Dt will be denoted by x If the direction of the outward normal is specified, then a circle may be uniquely defined by its radius and one point on its circumference. Consider then the circle passing

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

921

procedure is quite general and, in principle, can be applied to impose any contact angle. In the case of a dynamic contact angle, the value of the contact angle can be determined empirically using ‘‘Tanner’s law” (see [85,86]) uc ¼ F ðuÞ;

ð67Þ

where uc is the velocity of the contact point in the direction of the wall and u is the contact angle related to the surface and wall normals given by cosðuÞ ¼ n2  n3 :

ð68Þ

One simple approximation for F ðuÞ is given by F ðuÞ ¼ a0 þ a1 u

ð69Þ

for given constants a0 and a1 . ^2 is the contact point position at a previous Recall that x time, say t  Dt, the velocity of the contact point uc may be estimated using a first-order approximation uc ¼

Fig. 9. Method for the approximation of the contact angle.

through x1 and x2 with radius r. If the outward unit normal at x1 is n1 then this uniquely determines its centre which we shall denote by x0 . The two main identities follow immediately: x1  x0 ¼ rn1 ;

x2  x0 ¼ rn2 :

Therefore x1  x2 ¼ rðn1  n2 Þ:

ð63Þ

On the other hand, if b denotes the normal distance from x1 to the wall, we have that b ¼ ðx1  x2 Þ  n3 :

ð64Þ

^ 2 Þ  sw ðx2  x ; Dt

ð70Þ

where the direction of sw is determined by imposing sw  n2 > 0. The curvature, therefore, can be calculated by solving (65)–(69). A simple, but effective procedure is as follows. First compute b directly from (65). Use the old value of uc as an initial guess and compute u from (67) with (69). Use (65) and (68) to determine r ¼ 1=j and n2 . Eq. (66) is now checked to see if it is satisfied. If not, then a bisection (or other root finding) algorithm may be used until uc has been found to sufficient precision. Eq. (70) may then be used to calculate x2 . For low Reynolds number simulation, when Dt is small due to viscous stability restrictions, correction steps are not required.

Hence, on substituting (63) into (64), we obtain rðn1  n2 Þ  n3 ¼ b

ð65Þ

11.1. Data structure

and then j¼

1 ðn1  n2 Þ  n3 ¼ : r b

11. Data structure and system environment

ð66Þ

We do not explicitly impose a slip boundary condition at the contact line. For B-cells adjacent to a S-cell we apply standard no-slip boundary conditions, and that the computed capillary pressure is applied at the S-cell. This is due to the fact that the slip region is interpreted as restricted to a microscopic region in the proximity of the contact point, while in most of the B-cell the no-slip conditions will still prevail. Therefore, the position of the marker particles, which should be updated using the contact point velocity, will not be correct close to the wall. Hence, for the analysis of the results, the position of the surface given by the tracking points at the cell closest to the wall are only approximate and may be disregarded, while the effective extension of the free surface and the contact point are given by the approximation described in this Section. The above

Freeflow-3D uses two classes of data: direct and indirect data. The first class contains the data concerning the domain, velocities, pressure, cells, parameters used by the simulator (Simflow-3D) and the representation of geometric objects of the model. The direct data set is subdivided into two sub-classes, namely, static and dynamic data. The static data comprise the data defining the domain, the discretization, the non-dimensionalization parameters and the time-independent fluid properties (e.g. viscosity). The dynamic data comprise data concerning velocities, pressure, cell type and the representation of geometric objects. Velocity and pressure are stored as matrices. The cells are represented by a matrix and, additionally, those containing fluid, or on the rigid boundary, are also stored in a tree-like data structure. The reason for having some of the cells represented twice is because access to specific groups of cells, namely F, S, B and I cells, is required

922

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

throughout the code and in this case tree storage is much more efficient. On the other hand, when access to all cells or to adjacent cells is required a matrix representation fares better. Moreover, further information can be built into the tree representation, improving the performance of the code. For instance, for a S cell the normal vector needs to be stored, whereas for a boundary cell the intersections of the grid lines with the container surface are required. The geometric objects have their surfaces approximated by piecewise linear surfaces and so are represented by a half-edge data structure which is a B-Rep (Boundary Representation) [58] structure. This structure is used to represent faces, edges, vertices and their topological relationships. It was originally designed to model rigid solids from a knowledge of their surfaces. The insertion and deletion of data into the half-edge data structure are made by the Euler operators which also guarantee the consistency of the represented object. Among other types of representations of solids (e.g. Constructive Solid Geometry, Spatial Occupancy Enumeration, Cell Decomposition [58]), the B-Rep representation was chosen because of its conciseness, its explicit representation of the boundary surface and its capability to efficiently perform local modifications to the object’s surface. The data structure employed to represent the class of indirect data was constructed with three main points in mind: the representation of geometric objects by their surfaces, the performance of the algorithms involved in the simulation, and the minimization of interdependency of the data. The indirect data consists of two types, called Containers, Inflows and Outflows, representing containers, inflows and outflows, respectively. The Container data structure is composed of geometric data (B-Rep), boundary condition types, information stored in a tree structure about the cells that define the container and specific attributes of the container. In addition, the container can move by prescribing trajectories such as harmonic motion, rotational motion, piecewise linear trajectories, splines trajectories and movements defined by external forces. These types of movements are stored in the Container data structure. The Inflow data structure is composed of geometric data (B-Rep), inflow properties, information about the container and the fluid related to the inflow, and a tree structure storing information about the cells defining the inflow and specific attributes concerning the inflow. Like the containers, the inflows can move by prescribing the trajectories of harmonic motion, piecewise linear trajectories and splines trajectories. The information about the inflow movement is stored in the Inflow data structure. The Outflow data structure is also composed of geometric data (B-rep) and outflow properties such as its plane of definition and orientation. The data structures described above are used by Modflow-3D to build the initial configuration for the flow that will then be evolved by Simflow-3D.

11.2. The module Modflow-3D Modflow-3D is responsible for data initialization of the simulation module. It consists of an user-friendly graphic interface for data input related to the flow, such as the domain, mesh, cells, reference values, flow type, parameters, error tolerance used by the simulation module and geometric objects making up the model. Modflow-3D also permits the visualization of the entered data. The technique used for viewing objects is rendering which makes it easier for the user to interpret what he sees. In addition to the usual geometrical operations like rotation, translation and zoom, the user can choose the type of projection (parallel, perspective), the viewer’s distance, the rendering type (e.g. Wire-frame, Flat, Gouraud, Phong) and ambiental and directional light intensity. Another useful feature of this module is the possibility of altering object visualization properties like color, visibility, transparency and light reflection properties. The geometric objects are represented using the halfedge B-Rep data structure, and this representation is achieved using a well known modeling technique: threedimensional primitive instancing. The technique of primitive instancing permits the choice of objects from a menu of pre-defined objects by simply defining their parameters. Currently, there are several choices, such as open rectangular (trapezoidal) boxes, open cylindrical/conical boxes, solid sphere, solid cylinder, pipe, solid square torus and so on. The remaining data in the structures Container, Inflow, Outflow and Fluid are directly obtainable from a pulldown menu. For instance, the structure Inflow also contains information on: inflow plane and inflow orientation, inflow velocity, boundary condition type, parameters for particle insertion (number of particles in each direction in each cell), a tree storing the cell set defining the inflow discretization and specific attributes of the inflow. The inflows can be of type Linear or Parabolic. The inflow type Linear is defined by having a constant velocity profile while the Parabolic inflow has a parabolic velocity profile. Like the structure Inflow, the structure Outflow contains information about outflow plane and outflow orientation. The specific attributes are data such as: visualization parameters, and data related to the nozzle definition and movement. 11.3. The module Simflow-3D The module Simflow-3D is the central part of the simulation environment described in this article because it implements the governing equations and boundary conditions. One of the main challenges in generalizing the ideas of the two-dimensional code, GENSMAC [88], to three dimensions was in dealing with the free surface. This was because in two dimensions virtual particles were used to represent the fluid, and this technique could not have been carried over to three dimensions due to the very large number of particles that would have been needed to represent

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

the fluid to photographic precision. We resolved this problem by introducing a new procedure whereby marker particles were only employed on the surface of the fluid. This brought huge savings in storage and computing time making the extension of GENSMAC feasible for solving full three-dimensional problems efficiently. 11.4. The module Visflow-3D Visflow-3D is responsible for the presentation of the output produced by Modflow-3D and Simflow-3D. This output comprises data from the flow model, object representation and properties of the flow saved at pre-defined instants of time. Time-independent data are stored in a single file, while those which are time-dependent are stored in a sequence of files containing information about the geometry of objects, the flow properties, cell configuration, time and cycle number. Visflow-3D allows the user to view the geometry (containers, inflows and outflows), the fluid flow itself and also the flow properties (velocity, pressure, viscosity, shear rate and stress tensor components). Like Modflow-3D, Visflow-3D uses rendering techniques for object presentation in order to make user-understanding of the pictures more immediate. In addition to the viewing facilities provided by Modflow-3D, the user can view cuts of objects by planes parallel to one of the main axes or planes parallel to a user defined direction. Flow properties can be viewed in three dimensions using rendering techniques by considering the flow property as a texture, or they can be viewed as contour lines. Both the velocity and the stress tensor can be visualized through its components; for instance, the velocity can be understood by exhibiting its components u, v and w while the stress tensor can be interpreted by viewing its components S xx ; S xy ; . . . ; S zz . 11.5. Communication between the modules Communication between modules is carried out through files. Modflow-3D outputs six files: three of them store information on the half-edge structure of containers, inflows and fluids. Two files store information on the indirect structure of containers and inflows, and the sixth file stores all the remaining data such as viscosity, time-step factors, time-step size, scaling parameters, flow properties, etc. The module Simflow-3D has, as input, six files all with the same structure as the output files of the module Modflow-3D. These files may have been created either by Modflow-3D or by Simflow-3D. The output files of Simflow-3D have four different purposes: normal exit – creates files for continuing the simulation if desired; forced exit – files created if an error occurs, this aids debugging; pre-defined temporary saves – for the event of a system crash; and finally, visualization files – for viewing the results. In the first three cases the files have the same structure as those of Modflow-3D; in the last they have the structure of Visflow-3D. The input for the module Visflow-3D is created in

923

three files containing information on the half-edge structure of containers, inflows and fluids. These files are used for visualizing the objects. An additional five files, containing the values of the velocity and pressure fields u, v, w, p and the configuration of the cells, are created for visualizing the flow properties. A further file containing time-independent data (called the manager) is created at the beginning of the simulation and is used throughout. 12. Applications We conclude this paper by presenting numerical simulations of three-dimensional free surface flows. 12.1. Jet buckling When a cylindrical jet flows onto a rigid surface a phenomenon known as jet buckling can occur if the Reynolds number is smaller than a prescribed value. This flow has attracted the attention of a number of researchers and has been studied both experimentally and numerically by Cruickshank and Munson [15] and Cruickshank [16]. They presented results, both experimental and theoretical, for Newtonian jets. From their study, they obtained estimates for when jet buckling occurs. These are based on the Reynolds number and the aspect ratio H =D, where H is the height of the inlet to the rigid plate and D is the jet diameter. For an axisymmetric jet, they found that if the conditions Re < 1:2 and H =D > 7:0 were satisfied then the jet would buckle. In order to simulate this problem, we considered a cylindrical thin jet of diameter D ¼ 6 mm issuing from an inlet situated at a height H ¼ 5 cm above the plate (so that H =D 8:33) flowing onto a rigid plate. A uniform input velocity of U ¼ 1 m s1 was set at the inlet, and the fluid viscosity m ¼ 0:012 m2 s1 was used. Surface tension was neglected. The scaling parameters were U ; D; m so that Re ¼ 0:5 and 1=F 2r ¼ 0:242611. A mesh size of Dx ¼ Dy ¼ Dz ¼ 0:001 mm (85  85  155) was employed. Fig. 10 displays the fluid flow configuration at different times. Fig. 11 also displays the equivalent fluid flow configuration for a planar jet. 12.2. Splashing drop A spherical drop of fluid of diameter D ¼ 10 mm is given an initial velocity of U ¼ 1 m s1 and released from a height of 4 cm above a square container containing a quiescent fluid. The value of the viscosity was m ¼ 106 m2 s1 so that the Reynolds number Re ¼ UD=m ¼ 1; 000. The mesh used in this example was 100  100  100 cells (Dx ¼ Dy ¼ Dz ¼ 1:0 mm). The initial indentation may be observed with the formation of a traveling wave. High pressure beneath this indentation causes fluid to travel upwards (in the form of a splash). As the ‘‘splash” reaches its peak we can observe the surface waves being reflected from the sides of the container (see Fig. 12). In this

924

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

Fig. 12. Three-dimensional numerical simulation the splashing drop at different times. Fig. 10. Three-dimensional numerical simulation of cylindrical jet buckling at different times.

Fig. 13. Three-dimensional numerical simulation of a circular hydraulic jump at different times.

Fig. 11. Three-dimensional numerical simulation of planar jet buckling at different times.

problem surface tension is insignificant (see, for example, [51]) and was neglected. 12.3. Hydraulic jump The Freeflow3D code has been applied to simulate a circular hydraulic jump. We considered a cylindrical jet of a viscous fluid flowing rapidly onto a rigid horizontal surface. The non-dimensional parameters were: Re ¼ 400, Fr ¼ 1:785. The scaling parameters were: inflow velocity

U ¼ 0:5 m s1 , inflow diameter D ¼ 8:0 mm and kinematic viscosity m ¼ 105 m2 s1 . Gravity was considered to be acting in the z-direction with g ¼ 9:81 m s2 . The mesh used in this example was 100  100  100 cells (Dx ¼ Dy ¼ Dz ¼ 1:0 mm). The Freeflow3D code solved this problem with the above input data. Fig. 13 exhibits snapshots of this run at selected times. Surface tension has been neglected. 12.4. Rising bubbles This example displays the simulation of three phase fluid flow of two bubbles with different densities, viscosities and diameters, rising in a third phase with a greater density. The bubbles rise due to buoyancy forces. The non-dimensional parameters were: Re ¼ 55:75, Fr ¼ 1, xe ¼ 4:6. For 3 the continuous phase, qf ¼ 880 kg=m and lf ¼ 0:0125 N 2 3 s=m . For the larger bubble, D1 ¼ 4 mm, q1 ¼ 88 kg=m 2 and l1 ¼ 0:00125 N s=m , with Eo1 ¼ 4:6. Finally, for the smaller bubble, D2 ¼ 2 mm, q2 ¼ 176 kg=m3 and l2 ¼

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

925

0:0025 N s=m2 , with Eo2 ¼ 1:15. The surface tension coefficient for all interfaces was chosen to be rI ¼ 0:03 N=m. The lengths and velocities pffiffiffiffiffiffiffiffiwere non-dimensionalized using D1 ¼ 4 mm and U ¼ gD1 , respectively. The grid used in this example was 32  32  64 cells. Fig. 14 displays the time evolution of the flow. It can be seen that the larger bubble undergoes a larger deformation than the smaller bubble, due to its larger Eo¨tvo¨s number (Eo ¼ qgD2 r1 ). There is a small deflection in the trajectory of the smaller bubble as a result of the interaction with the larger bubble as it passes by. Surface tension plays an important role.

zle. The following input data were employed: The nondimensional parameters were: Re ¼ 7:5, Fr ¼ 2:9146. The scaling parameters were: velocity reference value U ¼ 1:0 m s1 (fluid velocity at the nozzle), length reference value L ¼ 12:0 mm (jet width in the x-direction), and kinematic viscosity 1:63 m2 s1 . Gravity was assumed to be acting in the z-direction with g ¼ 9:81 m s2 . The mesh used in this example was 80  120  70 cells (Dx ¼ Dy ¼ Dz ¼ 1:0 mm). The Freeflow3D code solved this problem with the above input data. Fig. 15 displays a series of snapshots taken during this run.

12.5. Container filling

13. Viscoelastic flows

We present a calculation which simulates the filling behavior of a trapezoidal container with a rectangular noz-

13.1. Oldroyd-B Complex fluids such as polymer melts and solutions, suspensions and micelle solutions are encountered in an variety of industries and applications critical to modern day economies, and a fundamental understanding of their behavior is becoming increasingly important. There are many models of viscoelasticity; the reader is referred to the recent book by Owens and Phillips [68]. In this Section we focus attention on just one, the Oldroyd-B model. This can be characterized by the equations of mass and momentum conservation (4) and (5) and its constitutive equation   5 5 s þ k1 s ¼ 2l D þ k2 D ; ð71Þ where k1 and k2 are the relaxation and retardation time constants. 13.2. Elastic jet

Fig. 14. Two bubbles rising in a continuous phase, at different times.

Fig. 15. Simulation of the filling process of a trapezoidal container at different times.

To illustrate the effect of viscoelasticity on the buckling phenomenon we applied Freeflow3D to simulate thin jets flowing onto a rigid plate. We considered a rigid plate of dimensions 6 cm  6 cm and an axisymmetric nozzle situated 12 cm above the rigid plate. The nozzle radius was R ¼ 3 mm and the velocity of the jet issuing from the nozzle was 1 m s1. Two simulations were performed: one run with a Newtonian jet and another run using a viscoelastic jet modeled by the Oldroyd-B constitutive equation. The fluid properties were q ¼ 1; 000 kg=m3 , l0 ¼ 0:004615 Pa s, k1 ¼ 0:006 s, k2 ¼ 0:0006 s. Therefore we had Re ¼ qUR=l0 ¼ 1:3 and We ¼ k1 U =R ¼ 1. A mesh size of 60  120 cells was employed in both simulations. We anticipated that the Newtonian jet will not buckle as the restriction Re < 1:2 is not satisfied (see Section 12.1). The results are displayed in Fig. 16 showing the jet flowing onto a rigid plate at different times. The Newtonian jet flows radially without any sign of the buckling instability while the jet modeled by the Oldroyd-B constitutive equation does undergo slight buckling. The reason why the Oldroyd-B jet buckles may be because of the extensional viscosity which increased after the jet impinges upon the rigid plate.

926

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

Fig. 16. Numerical simulation of jet impinging onto a rigid plate at different times. (a) Newtonian jet; (b) viscoelastic jet modeled by the Oldroyd-B model.

Fig. 17. Fluid flow visualization at different times of the simulation of the extrudate swell for increasing values of Weeffect .

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

13.3. Die-Swell We also applied the Freeflow3D code to the simulation of the extrudate swell of an Oldroyd-B fluid. The timedependent flow of an axisymmetric jet flowing inside a tube and then extruded in air was considered. The no-slip boundary condition is imposed on the tube walls, while fully developed flow is assumed at the tube entrance. On the fluid free surface, the full stress conditions are applied. To simulate this problem the following input data were employed: tube radius R ¼ 1 cm, tube length L ¼ 3R, Dx ¼ Dy ¼ Dz ¼ 0:1 cm (20  100  100 cells). The fluid was characterized by m0 ¼ 0:010 m2 s1 , k1 ¼ 0:01 s. The scaling parameters were R, U ¼ 1, m0 and k1 , giving Re ¼ UR=m0 ¼ 1 and We ¼ k1 U =R ¼ 1. To demonstrate that the code can deal with the Oldroyd-B model, we used these input data and performed three simulations. In the first simulation the value of k2 ¼ 0:9k1 was used and in the second simulation we used k2 ¼ 0:70k1 while in the third simulation we chose k2 ¼ 0:5k1 . Recall that the effective Weissenberg number for the Oldroyd-B model is given by (see [105]) Weeffect ¼ ð1  kk21 ÞWe. Thus, in these simulations we used Weeffect ¼ 0:1; 0:30; 0:5, respectively. Fig. 17 shows the jet flowing inside the tube and then being extruded into the air. For the time t ¼ 0:04s (first column) the jet is just leaving the tube and the differences between the three simulations are not noticeable, at least to the human eye. However, at later times the jet is extruded into the air and the differences between the three simulations become clearer. This is particularly true for the case of t ¼ 0:12 s (last column), where we can observe that the results for the case of Weeffect ¼ 0:5 present a much larger swelling than that of the other two simulations. Indeed, at the time t ¼ 0:12 s, the maximum swelling ratios (S r ¼ Rmax =R) for the three simulations were 28% for Weeffect ¼ 0:1, 46% for Weeffect ¼ 0:30 and 59% for Weeffect ¼ 0:5.

927

one such problem is that of bubbles rising in a beer glass – this could not be properly achieved with a billion computational cells (unless, of course, the bubbles were treated as a continuum). We have seem that in order to compute higher Reynolds number flows it was found necessary to employ high order monotone upwinding for the convective terms. In the future it will probably be necessary to sacrifice the essential simplicity of the MAC approach and introduce adaptivity and higher order methods or indeed other numerical methods such as finite elements or spectral collocation. The ultimate goal is virtual reality computational fluid dynamics (VRCFD) whereby the output is, ideally, in the form of a dynamic hologram (preferably in real time) which is essentially indistinguishable from the ‘‘real thing”. In Section 3 we presented a brief description of the volume of fluid approach, the level set method and several different front tracking techniques for solving free surface flow problems. Many codes have been written describing the many variants (and hybrid variants) of the above. These codes are necessarily complex and it is becoming increasingly difficult to compare one with the other, and thereby come to some kind of conclusion as to which is the best. Indeed the word ‘best’ may have little meaning here, as one code may be excellent at solving bubbly flow while another may be particularly effective when employed to compute viscoelastic flows. Nonetheless, it behoves the CFD community, in particular that large subset interested in free surface flows, to attempt a systematic interrogation of all existing codes (and their underlying methods). One way to do this is for the community to agree on a set of test problems that all good codes should be able to solve. Most of us have experienced using a commercial code that did not live up to the vendor’s claims: such codes, if subjected to the test set, could be pinpointed and discounted. To initiate this discussion the authors would like to conclude by proposing six test problems that all ‘good’ codes should be able to solve:

14. Future developments One of the attractive features of the MAC method is its simplicity: a fixed Eulerian grid is used – it does not ‘suffer’ from the complexity of adaptivity; finite difference methods are employed and these, for the most part, are at most second order; robustness is achieved through the use of small, often very small time steps especially in the case of low Reynolds number flows; the conjugate gradient solver tends to converge rapidly without the need for a preconditioner. Today, this ‘toy’ method of the 1960s has already become a useful scientific and technological tool, able to compute quite large three-dimensional problems on grids of the size of 1,000,000 computational cells using a modest computer and to display the results as a realistic video sequence. Tomorrow’s computers should make possible computations over 1,000,000,000 computational cells. However, greater computing power itself may not be sufficient. One can think of many problems where this is true:

1. 3D jet buckling (Newtonian and viscoelastic). 2. 3D splashing drop (Newtonian and viscoelastic): tungsten ball dropped from different heights – various configurations showing satellite drops from the resultant splash. 3. 3D extrudate swelling (Newtonian and viscoelastic). 4. 3D bursting bubble: a bubble bursting at the surface of a fluid giving rise to a jet or a splash from its base. 5. A large number of bubbles (>100, say) of different sizes, densities and viscosities rising (and interacting) in a fluid. 6. Rod climbing of a viscoelastic fluid at moderate Weissenberg number. We suspect that not many existing codes could solve these six problems. The MAC method, as described in this article, cannot as yet, but we can foresee it being able to do so in the relatively near future.

928

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

Acknowledgements We acknowledgment support from the FAPESP (Fundacßa˜o de Amparo a` Pesquisa do Estado de Sa˜o Paulo), Contract No. 04/16064-9. The second author has been supported by CNPq (Conselho Nacional de Desenvolvimento Cientı´ fico e Tecnolo´gico) under Grants Nos. 302048/ 2004-9, 504622/2003-0 and CAPES/GRICES Grant No. 136/5. References [1] Amsden A, Harlow F. The SMAC method: a numerical technique for calculating incompressible fluid flows. Technical Report LA4370, Los Alamos National Laboratory, 1970. [2] Aulisa E, Manservisi S, Scardovelli R, Zaleski S. A geometrical area-preserving volume-of-fluid advection method. J Comput Phys 2003;192:355–64. [3] Aulisa E, Manservisi S, Scardovelli R. A surface marker algorithm coupled to an area-preserving marker redistribution method for threedimensional interface tracking. J Comput Phys 2004;197:555–84. [4] Baker GR, Moore DW. The rise and distortion of a twodimensional gas bubble in an inviscid liquid. Phys Fluids A 1989;9:1451–9. [5] Bird RB, Armstrong RC, Hassager O. Dynamics of polymeric liquids. Fluid mechanics, 2nd ed., vol. I. New York: John Wiley & Sons; 1967. [6] Boulton-Stone JM, Blake JR. Gas bubbles bursting at a free surface. J Fluid Mech 1993;254:437–66. [7] Brackbill JU, Ruppel HM. FLIP: a method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions. J Comput Phys 1986;65:314–43. [8] Brackbill JU. The ringing instability in particle-in-cell calculations of low-speed flow. J Comput Phys 1988;75:469–92. [9] Brackbill JU. FLIP MHD: a particle-in-cell method for magnetohydrodynamics. J Comput Phys 1991;96:163–92. [10] Butler TD, O’Rourke PJ. A numerical method for two-dimensional unsteady reacting flows. In: Sixteenth international symposium on combustion. 1976. p. 1503–15. [11] Castelo A, Tome´ MF, Cesar CNL, McKee S, Cuminato JA. Freeflow-3D and integrated simulation system for three-dimensional free surface flows. J Comput Visual Sci 2000;2:199–210. [12] Chenn I-L, Glimm J, McBryan O, Plohr B, Yaniv S. Front tracking for gas dynamics. J Comput Phys 1986;62:83–110. [13] Craik AD, Latham RC, Fawkes MJ, Gribbon WF. The circular hydraulic jump. J Fluid Mech 1981;112:347–62. [14] Cross MM. Rheology of non-Newtonian fluids: a new flow equation for pseudo-plastic systems. J Coll Sci 1965;20:417–37. [15] Cruickshank JO, Munson BR. Viscous-fluid buckling of plane and axisymmetric jets. J Fluid Mech 1981;113:221–39. [16] Cruickshank JO. Low-Reynolds-number instabilities in stagnating jet flows. J Fluid Mech 1988;193:111–27. [17] Darip P, Glimm J, Lindquist B, Maesumi M, McBryan O. On the simulation of heterogeneous petroleum reservoirs. In: Wheeler M, editor. Numerical simulation in oil recovery. New York: Springer Verlag; 1988. [18] Deville MO. Numerical experiments on the MAC code for slow flow. J Comput Phys 1974;15:362–74. [19] Enright D, Fedkiw R, Ferziger J, Mitchell I. A hybrid particle level set method for improved interface capturing. J Comput Phys 2002;183:83–116. [20] Ervin EA, Tryggvason G. The rise of bubbles in a vertical shear flow. J Fluid Eng Trans ASME 1997;119:443–9. [21] Evans MW, Harlow FH. The particle-in-cell method for hydrodynamic calculations. Los Alamos National Laboratory Report LA2139, 1957.

[22] Fauci LJ, Peskin CS. A computational model of aquatic animal locomotion. J Comput Phys 1988;77:85–108. [23] Ferreira VG, Tome´ MF, Mangiavacchi N, Castelo A, Cuminato JA, Fortuna AO, et al. High order upwinding and the hydraulic jump. Int J Numer Methods Fluids 2002;39:549–83. [24] Ferreira VG, Mangiavacchi N, Tome´ MF, Castelo A, Cuminato JA, McKee S. Numerical simulation of turbulent free surface flow with two-equation k–e eddy-viscosity models. Int J Numer Methods Fluids 2004;44:347–75. [25] Fogelson AL, Peskin CS. A fast numerical method for solving the three-dimensional Stokes equations in the presence of gas suspended particles. J Comput Phys 1988;79:50–69. [26] Galaktionov OS, Anderson PD, Peters GWM, Van de Vosse FN. An adaptive front tracking technique for three-dimensional transient flows. Int J Numer Methods Fluids 2000;32:210–34. [27] Glimm J, Grove JW, Lindquist B, McBryan O, Tryggvason G. The bifurcation of tracked scalar waves. SIAM J Sci Comput 1988;9:61–79. [28] Glimm J, Grove JW, Li XL, Shyue KM, Zeng YN, Zhang Q. Threedimensional front tracking. SIAM J Sci Comput 1998;19:703–27. [29] Glimm J, Grove JW, Li XL, Tan DC. Robust computational algorithms for dynamic interface tracking in three dimensions. SIAM J Sci Comput 2000;21:2240–56. [30] Glimm J, Grove JW, Zhang Q. Interface tracking for axisymmetric flows. SIAM J Sci Comput 2002;24:208–36. [31] Glimm J, Li XL, Liu YJ, Xu ZL, Zhao N. Conservative front tracking with improved accuracy. SIAM J Sci Comput 2003;41:1926–47. [32] Golafshani M. A simple numerical technique for transient creep flows with free surfaces. Int J Numer Methods Fluids 1988;8:897–912. [33] Greaves D. A quadtree adaptive method for simulating fluid flows with moving interfaces. J Comput Phys 2004;194:35–56. [34] Han J, Tryggvason G. Secondary breakup of axisymmetric liquid drops. I. Acceleration by a constant body force. Phys Fluids 1999;11:3650–67. [35] Harlow F. hydrodynamic problems involving large fluid distortion. J Assoc Comput Mach 1957;4:137. [36] Harlow F, Welch JE. Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface. Phys Fluids 1965;8:2182–9. [37] Harlow F, Nakayama PI. Turbulence transport equations. Phys Fluids 1967;10:2323–32. [38] Harlow F, Nakayama PI. Transport of turbulence energy decay rate. Los Alamos National Laboratory LA-3854, 1968. [39] Harlow FH, Welch JE. A fast ICE solution procedure for flows with largely invariant compressiblity. J Comput Phys 1981;40:254–61. [40] Hirsch C. Numerical computation of internal and external flows, vols. 1 and 2. Wiley series in numerical methods in engineering, 1988. [41] Hirt CW, Shannon JP. Free-surface stress conditions for incompressible-flow calculations. J Comput Phys 1968;2:403–11. [42] Hirt CW, Nichols BD, Romero NC. SOLA–A numerical solution algorithm for transient fluid flows. Los Alamos National Laboratory LA-5852, 1972. [43] Hirt CW, Nichols BD. Volume of Fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39:201–25. [44] Hirt CW. Flow3D users manual. Flow Science Inc.; 1988. [45] Harvie DJE, Fletcher DF. A new volume of fluid advection algorithm: the defined donating region scheme. Int J Numer Methods Fluids 2001;35:151–72. [46] Jaeger M, Carin M. The front-tracking ALE method: application to a model of the freezing of cell suspensions. J Comput Phys 2002;179:704–35. [47] Jamet D, Torres D, Brackbill JU. On the theory and computation of surface tension: the elimination of parasitic currents through energy conservation in the second-gradient method. J Comput Phys 2002;182:262–76.

S. McKee et al. / Computers & Fluids 37 (2008) 907–930 [48] Kim MS, Lee WI. A new VOF-based numerical scheme for the simulation of fluid flow with free surface. Part I: New free surfacetracking algorithm and its verification. Int J Numer Methods Fluids 2003;42:765–90. [49] Kim MS, Park JS, Lee WI. A new VOF-based numerical scheme for the simulation of fluid flow with free surface. Part II: application to the cavity filling and sloshing problems. Int J Numer Methods Fluids 2003;42:791–812. [50] Kohno H, Tanahashi T. Numerical analysis of moving interfaces using a level set method coupled with adaptive mesh refinement. Int J Numer Methods Fluids 2004;45:921–44. [51] Lampe J, Dilalla R, Grimaldi J, Rothstein JP. Impact dynamics of drops on thin films of viscoelastic wormlike micelle solutions. J Non-Newton Fluid Mech 2005;125:11–23. [52] Lemos C. Higher-order schemes for free surface flows with arbitrary configurations. Int J Numer Methods Fluids 1996;23:545–66. [53] LeVeque RJ. Numerical methods for conservation laws. Lectures in mathematics. Zurich: ETH; 1992. [54] Li J, Renardy YY, Renardy M. Numerical simulation of breakup of a viscous drop in simple shear flow through a volume-of-fluid method. Phys Fluids 2000;12:269–82. [55] Lorstad D, Fuchs L. High-order surface tension VOF-model for 3D bubble flows with high density ratio. J Comput Phys 2004;200: 153–76. [56] Lorstad D, Francois M, Shyy W, Fuchs L. Assessment of volume of fluid and immersed boundary methods for droplet computations. Int J Numer Methods Fluids 2004;46:109–25. [57] Mangiavacchi N, Castelo A, Tome´ MF, Cuminato JA, Bambozzi de Oliveira ML, McKee S. An effective implementation of surface tension effects using the marker and cell methods for axisymmetric and planar flows. SIAM J Sci Comput 2005;26:1340–68. [58] Ma¨ntyla¨ M. An introduction to solid modeling. Rockville: Computer Science Press; 1988. [59] Markham G, Proctor MV. Modifications to the two-dimensional incompressible fluid flow code ZUNI to provide enhanced performance. C.E.G.B. report TPRD/L/0063/M82, 1983. [60] McQueen JF, Rutter P. Outline description of a recently implemented fluid flow code ZUNI. C.E.G.B. report LM/PHYS/258, 1981. [61] Miyata H, Nishimura S. Finite-difference simulation of nonlinear waves generated by ships of arbitrary three-dimensional configuration. J Comput Phys 1985;60:391–436. [62] Miyata H. Finite-difference simulation of breaking waves. J Comput Phys 1986;65:179–214. [63] Moretti G. Computation of flows with shocks. Annual Rev Fluid Mech 1987;19:313–37. [64] Nichols BD, Hirt CW. Improved free surface boundary conditions for numerical incompressible flow calculations. J Comput Phys 1971;8:434. [65] Nichols BD, Hirt CW, Hotchkins RS. SOLA-VOF: a solution algorithm for transient fluid flow with multiple free boundaries, Technical Report LA-8355, Los Alamos Scientific Laboratory, 1988. [66] Nourgaliev PR, Dinh TN, Theofanous TG. The characteristic-based matching (CBM) method for compressible flow with moving boundaries and interfaces. J Fluid Engng – Trans ASME 2004;126:586–604. [67] Oishi CM, Ferreira VG, Cuminato JA, Castelo A, Tome´ MF, Mangiavacchi N. Implementing implicit procedures in GENSMAC. TEMA (Tendeˆncias Mat Apl Comput) 2004;5:257–66. [68] Owens RG, Phillips TN. Computational rheology. Imperial College Press; 2002. [69] Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79:1–49. [70] Peskin CS. Numerical analysis of blood flow in the heart. J Comput Phys 1977;25:220–52.

929

[71] Popinet S, Zaleski S. A front-tracking algorithm for accurate representation of surface tension. Int J Numer Methods Fluids 1999;30:775–93. [72] Pracht WE. A Numerical method for calculating transient creeping flows. J Comput Phys 1971;7:46–60. [73] Quecedo M, Pastor M. Application of the level set method to the finite element solution of two-phases flows. Int J Numer Methods Eng 2001;50:645–63. [74] Ramshaw JD, Dukowicz JK. APACHE: a generalized-mesh Eulerian computer code for multicomponent chemically reactive fluid flow. Los Alamos National Laboratory Report LA-7427, 1979. [75] Rivard WC, Farmer OA, Butler TD. RICE: a computer program for multicomponent chemically reactive flows at all speeds. Los Alamos National Laboratory Report LA-5812, 1975. [76] Sethian JA. Theory, algorithms, and applications of level set methods for propagating interfaces. Acta Numerica. Cambridge, UK: Cambridge University Press; 1995. [77] Sethian JA. Level set methods: evolving interfaces in geometry. Fluid mechanics, computer vision and material sciences. Cambridge, UK: Cambridge University Press; 1996. [78] Shin S, Juric D. Modeling three-dimensional multiphase flow using a level contour reconstruction method for front tracking without connectivity. J Comput Phys 2002;180:427–70. [79] Sicilian JM, Hirt CW. An efficient computation scheme for tracking contaminant concentrations in fluid flows. J Comput Phys 1984;56:428–47. [80] Sousa FS, Mangiavacchi N, Castelo A, Nonato LG, Tome´ MF, Cuminato JA. Simulation of 3D free-surface flows with surface tension. In: Proceedings of 16th Brazilian congress of mechanical engineering (COBEM), Uberlaˆndia – MG, CD-ROM, 2001. [81] Sousa FS, Mangiavacchi N, Nonato LG, Castelo A, Tome´ MF, McKee S. A front-tracking/front-capturing method for the simulation of 3D multi-fluid flows with free surfaces. J Comput Phys 2004;198:469–99. [82] Sussman M, Smereka P. Axisymmetric free boundary problems. J Fluid Mech 1997;341:269–94. [83] Sussman M, Fatemi E. An efficient interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow. SIAM J Sci Comput 1999;20:1165–91. [84] Sussman M, Puckett EG. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible twophase flows. J Comput Phys 2000;162:301–37. [85] Tanner LH. The spreading of silicone oil on horizontal surfaces. J Phys D 1979;12:1473–84. [86] Tanner RI. Engineering Rheology. Cambridge, MA: Oxford University Press; 1992. [87] Theodorakakos A, Bergeles G. Simulation of sharp gas–liquid interface using VOF method and adaptive grid local refinement around the interface. Int J Numer Methods Fluids 2004;45:421–39. [88] Tome´ MF, McKee S. GENSMAC: a computational marker-andcell method for free surface flows in general domains. J Comput Phys 1994;110:171–86. [89] Tome´ MF, Duffy B, McKee S. A numerical technique for solving unsteady non-Newtonian free surface flows. J Non-Newton Fluid Mech 1996;62:9–34. [90] Tome´ MF, Castelo A, Cuminato JA, McKee S. GENSMAC3D: implementation of the Navier–Stokes equations and boundary conditions for 3D free surface flows, Universidade de Sa˜o Paulo, Departamento de Cieˆncia de Computacßa˜o e Estatı´stica, Notas do ICMC, No. 29, 1996. [91] Tome´ MF, McKee S. Numerical simulation of viscous flow: buckling of planar jets. Int J Numer Methods Fluids 1999;29:705–18. [92] Tome´ MF, McKee S, Barratt L, Jarvis DA, Patrick AJ. An experimental and numerical investigation of container filling: viscous liquids. Int J Numer Methods Fluids 1999;31:1333–53.

930

S. McKee et al. / Computers & Fluids 37 (2008) 907–930

[93] Tome´ MF, Castelo A, Cuminato JA, Murakami J, Minghim R, Oliveira CF, et al. Numerical simulation of axisymmetric free surface flows. J Comput Phys 2000;157:441–72. [94] Tome´ MF, Castelo A, Cuminato JA, McKee S. GENSMAC3D: a numerical method for solving unsteady three-dimensional free surfaces flows. Int J Numer Methods Fluids 2001;37:747–96. [95] Tome´ MF, Mangiavacchi N, Cuminato JA, Castelo A, McKee S. A numerical technique for solving unsteady viscoelastic free surface flows. J Non-Newton Fluid Mech 2002;106:61–106. [96] Tome´ MF, Grossi L, Castelo A, Cuminato JA, Mangiavacchi N, Ferreira VG, et al. A numerical method for solving three-dimensional generalized Newtonian free surface flows. J Non-Newton Fluid Mech 2004;123:85–103. [97] Tome´ MF, Castelo A, Federson F, Cuminato JA. A numerical method for solving the Oldroyd-B model for 3D free surface flows. In: Proceedings of XIV Congreso sobre Metodos Numericos Y Sys Aplicaciones – ENIEF 2004 – San Carlos de Bariloche, Argentina, 2004. [98] Torrey MD, Mjolsness RC, Stein LR. NASA-VOF3D: a three dimensional computer program for incompressible flows with free surfaces. Technical Report LA-11009-MS, Los Alamos National Laboratory, 1987. [99] Travis JR, Hirt CW, Rivard WC. Multidimensional effects in critical two-phase flow. Nucl Sci Eng 1978;68:338.

[100] Tryggvason G, Bunner B, Ebrat O, Tauber W. Computations of multiphase flows by a finite difference/front tracking method. I Multi-fluid flows. In: 29th Computational fluid dynamics. Lecture series 1998-03. Von Karman Institute for Fluid Dynamics; 1998. [101] Tsai TM, Miksis MJ. Dynamics of a drop in a constricted capillary tube. J Fluid Mech 1994;274:197–217. [102] Varonos A, Bergeles G. Development and assessment of a variableorder non-oscillatory scheme for convection term discretization. Int J Numer Methods Fluids 1998;26:1–16. [103] Wang JP, Borthwick AGL, Taylor RE. Finite-volume-type VOF method on dynamically adaptive quadtree grids. Int J Numer Methods Fluids 2004;45:485–508. [104] Welch JE, Harlow FH, Shannon JP, Daly BJ. The MAC method. Technical Report LA-3425, Los Alamos National Laboratory, 1965. [105] Yoo JY, Na Y. A numerical study of the planar contraction flow of a viscoelastic fluid using the simpler algorithm. J Non-Newton Fluid Mech 1991;30:89–106. [106] Yue WS, Lin CL, Patel VC. Numerical simulation of unsteady multidimensional free surface motions by level set method. Int J Numer Methods Fluids 2003;42:853–84. [107] Zhao Y, Tan HH, Zhang BL. A high-resolution characteristicsbased implicit dual time-stepping VOF method for free surface flow simulation on unstructured grids. J Comput Phys 2002;183: 233–73.