Numerical Simulation of Dispersed
Two-Phase Flow
PROEFSCHRIFT
ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,
op gezag van de Rector Magnificus prof. dr. ir. J.T. Fokkema, voorzitter van het College voor Promoties,
in het openbaar te verdedigen op dinsdag 18 september 2007 om 12:30 uur
door
Emil Rustam Ardeshir COYAJEE werktuigkundig ingenieur
Dit proefschrift is goedgekeurd door de promotoren: Prof. dr. ir. G. Ooms
Prof. dr. ir. P. Wesseling
Samenstelling promotiecommissie: Rector Magnificus, voorzitter
Prof. dr. ir. G. Ooms, Technische Universiteit Delft, promotor Prof. dr. ir. P. Wesseling, Technische Universiteit Delft, promotor
Prof. dr. ir. B.J. Boersma, Technische Universiteit Delft, toegevoegd promotor Prof. dr. A.E.P. Veldman , Rijksuniversiteit Groningen
Prof. dr. ir. J.A.M. Kuipers, Universiteit Twente
Prof. dr. R.F. Mudde, Technische Universiteit Delft Prof. R. Scardovelli, Universit´a di Bologna
The work presented in this thesis was supported financially by the
Netherlands Organisation for Scientific Research (NWO), project 635.000.004. Copyright c 2007 by E.R.A. Coyajee
Summary vii
Samenvatting ix
1 Introduction 1
1.1 Background and Motivation . . . 1
1.2 Computational Methods for Direct Numerical Simulation of In-terfacial Flow . . . 4
1.3 Objective . . . 7
1.4 Outline . . . 8
2 Governing Equations for Interfacial Two-Phase Flow 11 2.1 Introduction . . . 11
2.2 Fluid Motion . . . 11
2.3 Interface Description . . . 12
2.4 Dimensionless Parameters . . . 13
3 Single Marker Front-Capturing Method 17 3.1 Introduction . . . 17
3.2 Temporal Discretization of the Navier-Stokes Equations . . . 19
3.3 Spatial Discretization of the Navier-Stokes Equations . . . 21
3.4 Interface Advection . . . 24
3.4.1 CLSVOF Method . . . 25
3.4.2 MCLS Method . . . 28
3.4.3 Results of 2D Advection Test Cases . . . 31
3.4.4 Linear Advection of a Fluid Sphere . . . 40
3.4.5 Conclusion . . . 41
3.5 Interface Curvature . . . 42
3.6 Parallel Implementation . . . 45
3.6.1 Domain Decomposition Approach . . . 45
3.6.2 Conjugate Gradient Solver . . . 45
3.7 Results of the Single Marker Front-Capturing Method . . . 46
3.7.1 Stationary Droplet / Laplace Problem . . . 47
vi Contents
3.7.3 Buoyant Droplet in Creeping Flow . . . 50
3.7.4 Buoyant Droplet in Periodic Array at Moderate Reynolds Number . . . 52
3.8 Conclusion . . . 53
4 Multiple Marker Front-Capturing Method 63 4.1 Introduction . . . 63
4.2 Concept and Implementation . . . 65
4.3 Gravity-Driven Drop Impact on a Liquid-Liquid Interface . . . . 67
4.3.1 Experimental Set-up and Conditions . . . 68
4.3.2 Numerical Set-up . . . 69
4.3.3 Results . . . 70
4.4 Conclusion . . . 73
5 Level-Set Distance Function Reconstruction from Volume Frac-tions 87 5.1 Introduction . . . 88
5.2 Level-Set Interface Description . . . 90
5.3 Volume-of-Fluid Methodology . . . 91
5.4 Level-Set/Volume-of-Fluid Coupling in the MCLS Method . . . . 92
5.5 Review of Geometric Reconstructed Distance Function Methods 95 5.6 Reconstructing Distance Functions from Volume Fractions with the Fast Marching Method . . . 99
5.6.1 The Fast Marching Method . . . 100
5.6.2 Procedure for RDF from Volume Fractions with the Fast Marching Method . . . 101
5.6.3 Results . . . 104
5.7 Summary and Conclusion . . . 106
6 Conclusions and Recommendations for Future Work 111 A Volume Fraction Function 115 A.1 Forward Function . . . 115
A.2 Inverse Problem . . . 119
B The Eulerian-Lagrangian VOF advection scheme 121
Bibliography 129
Acknowledgments 131
A Front-Capturing Method for the Numerical Simulation of Dis-persed Two-Phase Flow - Emil Coyajee
A large variety of industrial processes involves a mixture of two immiscible liquids, e.g. emulsification in food industry and mixing in (petro-)chemical re-actors. Understanding the dynamics and interaction of bubbles and droplets in these processes is crucial for economically and ecologically optimized de-sign. Although experimental investigations provide essential insight into the phenomenology of such dispersed two-phase flows, obtaining quantitative data is difficult. Therefore, in the present work a computational method is developed for the detailed study of dispersed gas-liquid and liquid-liquid flows.
The computational method uses the approach of direct numerical simulation for interfacial two-phase flow. This approach is based on a direct mathematical description of the phases, considering fluid motion both inside the bubbles or drops and in the surrounding liquid. The boundary between the fluid phases is regarded as a moving, deformable interface. Each fluid phase has its own density and viscosity and the effect of surface tension is taken into account at the deformable interface.
The equations describing interfacial two-phase flow are solved with a front-capturing technique. The fluid motion is computed on a fixed grid, while a separate procedure is used to describe the geometry and the evolution of the interface. Moreover, jump conditions for the discontinuous variation of material properties and surface tension at the moving interface need to be incorporated into the spatial discretization of the equations for fluid motion on the fixed grid. In the development of the front-capturing method, specific attention has been directed to the following issues: The accurate representation and advection of the interface geometry, the accurate and robust implementation of interfacial jump conditions, a parallel implementation to allow large scale computations and the validation of the method. Furthermore, a fundamental issue in the modeling of two-phase flow concerns the collision of bubbles or drops. In reality, the collision of bubbles or drops may lead to coalescence. This depends not only on flow conditions, but also on physical phenomena on a molecular scale. Because the current mathematical model for the description of fluid motion is based on the continuum hypothesis, coalescence of bubbles or droplets is excluded from the numerical method.
viii Summary method is used, which has the advantage of the Volume-of-Fluid methodology for volume-conserving advection of the fluid phases and the Level-Set interface description for implementation of interfacial jump conditions. The properties of several recently developed combined Level-Set/Volume-of-Fluid methods are investigated in detail. For the implementation of interfacial jump conditions, also a hybrid approach is applied. At the interface between fluids, a continuous transition of viscosity is imposed, while a sharp implementation of the jump in the density and the pressure due to surface tension is achieved by means of the Ghost-Fluid methodology.
In the traditional Level-Set and Volume-of-Fluid methodology, reconnection of interfaces is automatically dealt with and coalescence of droplets occurs when-ever the distance between the interfaces is reduced to or below the mesh width. To prevent the numerical, artificial coalescence, a procedure has been developed which enables the representation of multiple interfaces in a computational cell within the front-capturing methodology.
A parallel implementation of the method has been realized to enable large scale computations for transient, three-dimensional studies. A domain decom-positioning approach is used for the parallelization.
Een ‘Front-Capturing’ Methode voor de Numerieke Simulatie van Gedispergeerde Tweefase Stroming - Emil Coyajee
Stromingen van twee niet-mengbare vloeistoffen spelen een belangrijke rol in een grote verscheidenheid van industri¨ele processen, bijvoorbeeld emulsifi-catie in de voedselindustrie en mengprocessen in (petro-)chemische reactoren. Begrip van de dynamica en de interactie van bellen en druppels in deze pro-cessen is essentieel om een economisch en ecologisch verantwoord ontwerp te kunnen realiseren. Alhoewel experimentele studies een goed inzicht geven in de fenomenologie van dit soort gedispergeerde twee-fase stromingen, is het moeilijk quantitatieve gegevens te verkrijgen. In het huidige werk is daarom een bereke-ningsmethode ontwikkeld die een zeer gedetailleerde studie van stromingen in twee-fase dispersies mogelijk maakt.
De berekeningsmethode is gebaseerd op een aanpak die bekend staat als directe numerieke simulatie van twee-fase stroming met een grensvlak. Deze aanpak is gebaseerd op een gedetailleerde wiskundige beschrijving van de fases, waarbij de vloeistofstroming zowel binnen de bellen en druppels als in de om-liggende vloeistof beschouwd wordt. De overgang van de ene naar de andere vloeistof fase wordt beschouwd als een vrij bewegend en vervormend grensvlak. Elke vloeistof fase heeft een eigen dichtheid en viscositeit en op het grensvlak wordt rekening gehouden met het effect van oppervlaktespanning.
De vergelijkingen die twee-fase stroming met een grensvlak beschrijven wor-den opgelost volgens de zogenaamde ‘front-capturing’ methode. De vloeistof-stroming van beide fases wordt berekend op een stationair grid, terwijl een aparte procedure wordt gebruikt om de geometrie en de evolutie van het grens-vlak tussen de fases te beschrijven. Bovendien moeten voorwaarden voor de discontinue variatie van de dichtheid, de viscositeit en de oppervlakte spanning aan het bewegende grensvlak opgenomen worden in de ruimtelijke discretisatie van de stromingsvergelijkingen op het stationaire grid.
In de ontwikkeling van de ‘front-capturing’ methode is speciale aandacht besteed aan de volgende zaken: De nauwkeurige representatie en advectie van de grensvlak geometrie, de nauwkeurige en robuuste implementatie van de voor-waarden aan het grensvlak, een parallelle implementatie van de methode om grootschalige berekeningen mogelijk te maken en de validatie van de methode. Bovendien is er een fundamentele kwestie in de modellering van twee-fase stro-mingen die de botsing van bellen of druppels betreft. In realiteit kan de
x Samenvatting ing van bellen of druppels ertoe leiden dat deze samensmelten, oftewel coa-lesceren. Dit hangt niet alleen af van de stromingscondities, maar ook van fysische processen op moleculaire schaal. Omdat in het huidige wiskundige model de vloeistoffen beschreven zijn op grond van de continuum hypothese, is coalescentie buiten beschouwing gelaten in de numerieke methode.
Voor de beschrijving van het grensvlak tussen de vloeistoffen is een gecom-bineerde ‘Level-Set/Volume-of-Fluid’ methode gebruikt. Dit heeft het voordeel dat de ‘Volume-of-Fluid’ methodologie gebruikt kan worden voor de volume be-houdende advectie van de vloeistof fases en de ‘Level-Set’ beschrijving van de grensvlak geometrie voor de implementatie van de eerdergenoemde voorwaar-den aan het grensvlak. De eigenschappen van verschillende recent ontwikkelde ‘Level-Set/Volume-of-Fluid’ methodes is uitgebreid onderzocht. Net als voor de beschrijving van het grensvlak zelf, wordt ook voor de implementatie van de voorwaarden aan het grensvlak een gecombineerde aanpak gebruikt. Ter plaatse van het grensvlak wordt een continue overgang van de viscositeit opgelegd, ter-wijl een scherpe overgang wordt gerealiseerd voor de sprong van de dichtheid en van de druk ten gevolge van oppervlakte spanning door middel van de ‘Ghost-Fluid’ methode.
In de traditionele ‘Level-Set’ en ‘Volume-of-Fluid’ methodes wordt recon-nectie van grensvlakken automatisch bewerkstelligd en coalescentie van drup-pels vindt plaats, zodra de afstand tussen de grensvlakken kleiner wordt dan de afstand tussen twee gridpunten. Om deze kunstmatige, numerieke coalescentie te verhinderen is een procedure ontwikkeld die de representatie van meerdere grensvlakken in een enkele grid-cel mogelijk maakt binnen de ‘front-capturing’ methodologie.
Introduction
1.1
Background and Motivation
Two-phase flow is a vast topic within the research area of fluid dynamics. This thesis specifically deals with two-phase flows in which a fluid, gas or liquid, is dispersed respectively as bubbles or drops in a continuous liquid. Such flows oc-cur in a large variety of industrial processes. Examples are emulsification in food industry, mixing in (petro-)chemical reactors and separation in oil production. Understanding the dynamics and interaction of bubbles and droplets in these processes is critical for optimization of design, i.e. maximizing the efficiency of the process while minimizing the costs and environmental impact. Although the behavior of single droplets is known to a large extent, much less is understood of the dynamics and interaction of multiple bubbles or droplets in a liquid. To improve our understanding for three kind of flows, research can be conducted in a number of ways.
On the one hand, experimental studies are indispensable for the investiga-tion of the behavior of dispersed liquid-liquid or gas-liquid flows. Experimental visualizations offer the possibility of getting insight into the phenomenology of the flow. For instance, the relation between flow conditions and flow regime can be studied and phenomena such as break-up and coalescence of droplets and bubbles can be observed. An example is the experimental study of an oil-water mixture in a turbulent Couette flow by Ravelet, Delfos and Westerweel (2007). A schematic of the set-up and a photograph of an instantaneous distribution of fluid phases in this experiment is shown in figure 1.1.
Such a visualization inside a mixture of two liquids is difficult to obtain, due to optical difficulties. When the refractive index of the two phases is not identical, the curved interface between the two fluids acts as a lens, providing a distorted image of the field of view. For the same reason, it can be difficult to extract detailed quantitative results from measurements. In relatively simple systems, e.g. the rise of a single droplet or the collision of two drops, accurate results can be obtained: By matching the refractive indices of the fluids and
2 Chapter 1 r = r r = 120 mm ωi= 2 Hz o ω r = 110 mm g h = 220 mm o i iωi oωo
Figure 1.1: Experimental set-up and snapshot of an oil-water dispersion in turbulent Couette flow (Ravelet et al. 2007).
adding tiny tracer particles to the liquid phases, the velocity field both inside the droplet and the surrounding liquid can be measured using Particle Image Velocimetry (PIV) (e.g. Mohamed-Kassim and Longmire 2003, Poesio, Poelma, Kaijen, Delfos and Westerweel 2007). In case of more complex flows, such as the turbulent Couette flow of an oil-water mixture in figure 1.1, it is still possible to configure the liquids so that refractive indices are matched. However, because of mixing induced by the turbulence, tracer particles are continuously moved around, also toward the interface between the two-phases. Here the particles tend to accumulate, making visualizations and PIV measurements difficult, if not impossible.
Because of the limitations in obtaining quantitative data from measurements and the costs involved with the construction and execution of experiments, com-puter simulations conducted with numerical models offer an important alterna-tive to: (1) Improve our understanding of dispersed fluid-liquid flows, (2) Assist in the modeling and design of two-phase processes in industry.
volume−averaged description Continuum approach using a
of phases Design industrial flows Design industrial flows/ Research tool
Research tool
Two−Fluid Method
Representation of bubbles/ droplets as discrete particles in the continuous liquid
Direct description of each individual phase
DNS of Interfacial Flow Discrete Bubble Method
smaller scales
larger geometry
Figure 1.2: Classification of computational techniques for dispersed gas-liquid/liquid-liquid flows. After Kuipers (2003).
suit different applications. A common classification of computational methods for dispersed gas-liquid or liquid-liquid flows is given by subdivision into the following categories: The Two-Fluid approach, the Discrete Bubble method and methods for the Direct Numerical Simulation (DNS) of interfacial two-phase flow (figure 1.2).
In the Two-Fluid and Discrete Bubble methods, a simplified description is used for the representation of bubbles or droplets in the continuous liquid. Therefore, these methods are relatively affordable in terms of computational effort and often applied for the modeling and design of large scale industrial flows (Portela and Oliemans 2006). Due to the simplified description however, information of the small flow scales is lost. This leads to the so-called ‘closure problem’: Additional terms are required in the equations for fluid motion, e.g. to represent the interaction between droplets or bubbles and the continuous liq-uid. The performance of the Two-Fluid and Discrete Bubble methods critically depends on the quality of the closure terms. Unfortunately the accurate defi-nition of these terms is a challenge on its own and the performance of closure models may vary with the complexity of a particular problem.
4 Chapter 1 number of drops or bubbles and the flow conditions in a simulation, e.g. the (droplet) Reynolds number, are limited.
In the present work, the development of a method for the DNS of dispersed gas-liquid and liquid-liquid flows will be considered. The approach is motivated by a desire to obtain detailed results for the behavior and interaction of multiple bubbles or droplets in a liquid, which are difficult to obtain from experiments or impossible to extract from simulations using Two-Fluid and Discrete Bubble models. Given the choice for such a demanding method, its development requires specific attention to enable not only accurate, but also efficient simulation of a large amount of droplets using parallel computing facilities.
1.2
Computational Methods for Direct
Numer-ical Simulation of Interfacial Flow
Efforts to compute free-surface and interfacial flows go back to the first develop-ments in computational fluid dynamics (e.g. Harlow and Welch 1965). However, solving the full three-dimensional Navier-Stokes equations with variable density and viscosity in the presence of a deformable interface with surface tension is a highly demanding task. Therefore, progress in the area of interfacial multiphase flow has been slow in comparison with the development of simulation tools for single-phase flows, which are now widely used both for research and engineering purposes. Until recently, simulations of three-dimensional interfacial flows at fi-nite Reynolds number, i.e. not concerning the asymptotic limits of the inviscid or creeping flow regime, were limited to very simple problems.
Compared to DNS of single-phase flow, additional challenges arise in the numerical modeling of two-phase flow with deformable interfaces:
• Representation of interface jump conditions. On a continuum scale, the material properties (density and viscosity) vary discontinuously at the interface and surface tension needs to be accounted for as a force concen-trated at the interface. Mathematically, the discontinuous viscosity and surface tension are incorporated into a set of jump conditions for tangen-tial and normal stresses at the interface.
• Representation and evolution of the interface geometry. Here, a notori-ous issue is the accurate estimation of interface normals and curvature for the implementation of jump conditions. Another issue is the mass-conserving advection of the fluid phases in case of incompressible flow. Finally, topology change and reconnection of interfaces is required when considering break-up and coalescence of droplets.
us-<0 =0 >0 φ φ φ
Figure 1.3: Front-tracking methodol-ogy: Discrete approximation of the in-terface with marker particles.
Figure 1.4: Level-Set methodology: Interface definition by iso-surface of a marker function, φ, defined at discrete locations on the fixed grid (crosses).
ing: (1) Deformable grids, locally adapted to the interface, or (2) fixed grids, using a separate procedure to describe the position of the interface.
The use of a deformable grid is appealing because of the straightforward def-inition of interface geometry by the boundaries of the discrete grid elements and the accurate representation of interface jump conditions. However, large defor-mations of the grid may occur in simulations with a large number of randomly moving, arbitrarily shaped drops. In this case, complete re-meshing of the com-putational domain may be necessary at regular intervals to maintain accurate representation of the interfaces. Therefore, the application of a deformable grid for the DNS of a two-phase dispersion can become computationally complex or expensive and will not be considered in the present work.
6 Chapter 1 0.02 0.01 0 1 1 1 0.5 0.3 0.7
(a) Analytical volume fractions for a circular interface on the fixed grid.
(b) Piecewise linear interface construc-tion from volume fracconstruc-tions.
Figure 1.5: Volume-of-Fluid methodology.
Implicit interface representations use a separate phase function, discretized on the fixed grid, indicating which phase is present at a given location. The interface is then implicitly defined by the location where the phase indicator function changes property. Therefore, this approach is also referred to as front-capturing. Two classical front-capturing approaches are the Volume-of-Fluid and the Level-Set technique, see Scardovelli and Zaleski (1999), Osher and Fed-kiw (2001) for reviews.
In the Level-Set methodology, the interface is represented by an iso-surface (usually the zero-level) of a smooth indicator function, commonly referred to as the Level-Set function (figure 1.4). Away from the interface, the Level-Set function is typically maintained a distance function to the interface. The main advantage of this representation is that the interface location, normal and cur-vature can be easily extracted from Level-Set grid values using finite difference approximations to partial differential equations. This is very convenient for the implementation of interface jump conditions. The main disadvantage of the Level-Set description is the difficulty of ensuring volume-conserving advection of fluid phases. Moreover, changes in interface topology are automatically dealt with in the Level-Set method. Although the Level-Set function itself can be dis-cretely conserved during advection, this does not necessarily imply conservation of the volume enclosed by the Level-Set iso-surface.
We conclude that the classical front-tracking, Volume-of-Fluid and Level-Set methodologies each have different advantages and disadvantages. Neither of them can be identified as ‘the ideal methodology’ for meeting the previ-ously defined challenges in the DNS of interfacial flow. As a result, new vari-ants and improvements are frequently published, often combining character-istics of several methods. Favorable combinations are either the Level-Set or front-tracking approach with the Volume-of-Fluid methodology. For these com-binations the advantages and disadvantages of the individual methods, con-cerning interface representation and volume conserving interface advection, are complementary. An example is for instance the mixed markers and Volume-of-Fluid method by Aulisa, Manservisi and Scardovelli (2003), combining the front-tracking and VOF methodology. A combined Level-Set and Volume-of-Fluid method has been pursued by various authors, e.g. the Coupled Level-Set/Volume-of-Fluid (CLSVOF) method by Sussman and Puckett (2000) and the recent Mass-Conserving Level-Set (MCLS) method by Van der Pijl (2005).
1.3
Objective
The objective of this work is the development of a front-capturing method for the direct numerical simulation of two-phase dispersions, i.e. multiple bubbles or droplets in liquid flow. A fixed grid approach will be used for the computation of the fluid motion, while the interface description will be based on a combined Level-Set/Volume-of-Fluid method. To ensure the accuracy and efficiency of the method, specific attention in its development will be given to the following: • Comparison and evaluation of different combined
Level-Set/Volume-of-Fluid methods for interface representation and evolution.
• Comparison and evaluation of different approaches for implementation of interface jump conditions.
• Parallelization of the method to enable large scale computations.
• Validation of the method against analytical solutions and results from experimental investigations.
Finally, a fundamental issue needs to be considered, concerning the colli-sion of bubbles or droplets. In reality, the collicolli-sion of different fluid bodies may lead to coalescence. This depends not only on flow conditions, but also on physical phenomena on a molecular scale. Because we will just consider the modeling of fluid motion based on the continuum hypothesis, we have decided to exclude coalescence of bubbles or droplets from the current numerical me-thod. Consequently, only dispersions of a fixed droplet size distribution can be investigated, but this appears to be a suitable starting point for investigating two-phase dispersions.
8 Chapter 1 <0 φ >0 φ >0 φ =0 φ φ=0 >0 φ <0 φ
(a) Before reconnection of interfaces. (b) After reconnection of interfaces. Figure 1.6: Numerical coalescence of approaching fluid bodies, represented by the Level-Set phase indicator function φ, which is defined at discrete locations on the fixed grid (crosses).
the distance between the interfaces is reduced to or below the mesh width. Figure 1.6 gives an example of such a situation for the Level-Set method. Obvi-ously, numerical coalescence can be avoided by extensive local grid refinement between different interfaces. However, local grid refinement increases compu-tational complexity and expenses considerably, especially in simulations with multiple, simultaneously colliding droplets. Therefore, a novel and economical approach will be proposed for the representation of colliding interfaces without coalescence within a combined Level-Set/Volume-of-Fluid method.
1.4
Outline
The content of thesis can be divided in two main parts, which can be read independently if desired.
their performance against several interface advection test-cases. After discussing the parallelization of the method with a domain decomposition approach, chap-ter 3 concludes with a validation of the complete front-capturing method. Its performance is evaluated for several simple problems, including some for which analytical solutions have been derived. In chapter 4, the single marker front-capturing method is extended to deal with multiple colliding, non-coalescing interfaces. The novel multiple marker front-capturing method will be validated by its application to the problem of gravity-driven drop impact on a liquid-liquid interface. Numerical results are directly compared to the results of an experimental study by Mohamed-Kassim and Longmire (2003).
The second part, which consists of chapter 5, concerns a broader study of combined Level-Set/Volume-of-Fluid methods, which focuses on the coupling procedure for the Level-Set interface and Volume-of-Fluid phase representations. This work is motivated by some of the less agreeable results that were obtained during the development of our front-capturing method. These results indicated that the coupling procedure is a crucial component in the performance of the combined Level-Set/Volume-of-Fluid methods, sometimes enhancing unphysical flow phenomena.
Governing Equations for
Interfacial Two-Phase Flow
In this chapter, the governing equations for interfacial two-phase flow are pre-sented. Subsequently, dimensionless parameters are derived by scaling the gov-erning equations. Finally, the relation between dimensionless parameters and shape regimes of individual bubbles and droplets is illustrated by means of the graphical correlation by Grace, Wairegi and Nguyen (1976).2.1
Introduction
In this work we consider the motion of two immiscible fluids, one fluid phase being dispersed in a surrounding continuous phase. Both fluid phases are allowed to have different density and viscosity, but these material properties are assumed to be constant within each phase. This is a simplification of reality, because fluid properties vary continuously over a thin layer of molecular dimension at the interface between the phases. On the continuum scale, this transition zone is regarded to be infinitesimally thin. Consequently, in the current mathematical model a discontinuous variation of density and viscosity from one bulk phase to the other is considered. Furthermore, the effect of surface tension is accounted for as a force concentrated at the interface.
Following the front-capturing methodology, two-phase flow is described by a set of equations for fluid and interface motion.
2.2
Fluid Motion
Both fluid phases are considered to be incompressible Newtonian fluids of con-stant temperature. In this case, the motion of the fluids is governed by the following equations (Batchelor 1974, Landau and Lifshitz 1987), describing
12 Chapter 2 Γ µ ρ, 2 µ ,ρ1 Phase ’1’ κn σ Phase ’2’ 1 2
Figure 2.1: Continuous fluid phase ‘1’ and dispersed fluid phase ‘2’ with cor-responding densities ρ1, ρ2 and viscosities µ1, µ2. Surface tension acts as a
force in the direction of the normal vector n at interface Γ and its magnitude is proportional to the surface tension coefficient σ and the local curvature κ. servation of mass:
∇ · u = 0, (2.1)
and conservation of momentum:
ρ (ut+ ∇ · (uu)) = −∇p + ∇ · µ ∇u + ∇uT + ρg, (2.2)
the latter equation also being referred to as the Navier-Stokes equation. Here p, u = (u, v, w)T and g respectively denote pressure, the velocity vector and the
gravitational acceleration. Subscript ‘t’ denotes differentiation with respect to time. The density ρ and viscosity µ are defined by ρ1, µ1in fluid phase ‘1’ and
by ρ2, µ2 in fluid phase ‘2’ (figure 2.1).
Since the fluids are considered immiscible, their mutual interface is a material property of the flow and the following conditions for continuity of fluid velocity and continuity of stresses at the interface apply (Batchelor 1974, Landau and Lifshitz 1987):
[u]Γ = 0, (2.3)
[pn + µ ∇u + ∇uT · n]Γ= σκn. (2.4)
Here [.]Γ denotes a jump across the interface Γ, n denotes the interface normal
vector, κ represents the magnitude of the interface curvature, and σ is the surface tension coefficient, which is assumed to be constant.
2.3
Interface Description
is denoted by φ. Thus, the mathematical description of the interface, denoted by Γ, is given by Γ(t) = {x | φ(x, t) = 0}. Away from the interface, φ, the so-called Level-Set function, is required to be a distance function to the interface such that φ < 0 in phase ‘1’ and φ > 0 in phase ‘2’. Consequently, the density and viscosity in (2.2) are defined by:
ρ = ρ1(1 − H(φ)) + ρ2H(φ), (2.5)
µ = µ1(1 − H(φ)) + µ2H(φ), (2.6)
where the Heaviside function H(φ) is defined by: H(φ) =
1 if φ > 0,
0 otherwise. (2.7)
Because the interface is a material property of the flow, its motion is described in terms of the Level-Set function by means of the following equation:
φt+ u · ∇φ = 0. (2.8)
As noted in chapter 1, numerical solutions of (2.8) do not conserve the volume enclosed by the Level-Set function, which is a requirement in the simulation of incompressible two-phase flow. Therefore, (2.8) is solved in combination with a volume-conserving approach for phase transport, the Volume-of-Fluid method. The Volume-of-Fluid methodology and the numerical procedure for interface advection are presented in section 3.4.
2.4
Dimensionless Parameters
For future use, scaling of the equations is introduced, leading to a number of relevant dimensionless parameters characteristic of two-phase flow. Dimension-less variables are obtained using characteristic length scale L, velocity scale U , density ρ1, and viscosity µ1:
x0 = x/L, u0
= u/U, t0 = t/(L/U ), ρ0
= ρ/ρ1, µ0 = µ/µ1, p0 = p/(ρ1U2).
Inserting these variables into equations (2.2) and (2.4) and dropping primes, we obtain the following dimensionless Navier-Stokes equation:
ρ (ut+ ∇ · (uu)) = −∇p +
1
Re∇ · µ ∇u + ∇u
T + ρ
F rz, (2.9) and condition for continuity of stresses at the interface:
14 Chapter 2 Here, z is a unit vector in the direction of the gravitational acceleration and the dimensionless density and viscosity are defined by:
ρ = 1 + (ζ − 1) H(φ), (2.11)
µ = 1 + (λ − 1) H(φ), (2.12)
where the density and viscosity ratios are respectively defined by ζ = ρ2/ρ1and
λ = µ2/µ1. The dimensionless parameters in (2.9) and (2.10) are the Reynolds
number, Weber number and Froude number, respectively defined by:
Re = ρ1U L/µ1, W e = ρ1LU2/σ, F r = U2/gL. (2.13)
The Reynolds, Weber and Froude number are measures for the relative impor-tance of inertia compared to respectively viscous, surface tension, and gravita-tional forces.
When a flow of two liquids is considered, the difference in densities may be relatively small compared to the densities of the individual phases. In this case a relevant definition of the Froude number is:
F r = ρ1 ∆ρ
U2
gL, (2.14)
where ∆ρ is defined by the magnitude of the difference between the densities of the phases. By taking the ratio of the Froude and the Weber number, the E¨otv¨os number is obtained:
Eo = ∆ρgL
2
σ , (2.15)
which is a measure for the relative importance of buoyancy/gravitational and surface tension forces.
A useful classification of the behavior of individual droplets and bubbles in terms of dimensionless parameters is presented in figure 2.2. This figure concerns the generalized graphical correlation by Grace et al. (1976) for bubbles and droplets rising or settling in an unbounded domain, which is based on just three dimensionless parameters. These are the droplet Reynolds number, Red,
the E¨otv¨os number, Eo, and the Morton number, M , defined by: Red=µcUdD ρc , Eo = ∆ρgD 2 σ , M = gµ4 c∆ρ ρ2 cσ3 . (2.16)
Here, Ud denotes the droplet terminal velocity, D denotes the sphere-equivalent
Single Marker
Front-Capturing Method
In this chapter, the computational approach is presented in case that all in-terfaces between different fluid phases are described by a single marker func-tion. First, the time integration procedure and the spatial discretization of the equations for fluid motion are discussed. For the spatial discretization, the im-plementation of surface tension is considered both with a sharp or continuum approach, i.e. the Ghost Fluid or the Continuous Surface Force method. Next, a concise description is given of two methods which are considered for the advec-tion of the interface: The Coupled Level-Set VoluOf-Fluid (CLSVOF) me-thod (Sussman and Puckett 2000) and the Mass Conserving Level-Set (MCLS) method (Van der Pijl 2005). The selection of the CLSVOF method is motivated by directly comparing the performance of the CLSVOF and MCLS methods for several test cases. Subsequently, the interface curvature representation is dis-cussed, followed by the approach for parallel implementation of the code.
Finally, results of the single marker front-capturing method are presented for a number of benchmark test cases for two-phase flows with surface tension, e.g. the stationary bubble (Laplace) test, a capillary wave between two superposed fluids, a buoyant droplet in creeping flow conditions and a buoyant droplet in periodic array at moderate Reynolds number.
3.1
Introduction
In the front-capturing methodology, numerical simulation of two-phase flow is performed as follows: The equations of fluid motion (2.1), (2.2) are solved on a fixed computational grid, while the material properties of the fluids and the jump condition for the stress are ‘captured’, i.e. accounted for in the spatial discretization of the Navier-Stokes equation, by tracking the location of the in-terface as it moves with the flow. As described in chapter 1, two main challenges of the numerical computation of interfacial two-phase flow are the description
18 Chapter 3 and evolution of the interface and the representation of discontinuous material properties and jump conditions for stress at the interface.
To deal with the discontinuous density and viscosity (2.5), (2.6) and the jump condition for the stress (2.4) in the spatial discretization of the Navier-Stokes equation, different approaches can be used: (1) Enforcing continuity of the solution by regularization (i.e. smoothing) of material properties and jump conditions; (2) Discontinuous representation of material properties and direct incorporation of jump conditions in the evaluation of spatial derivatives.
The first approach is also known as the Continuous Surface Force method-ology (CSF) (Brackbill, Kothe and Zemach 1992) or Continuous Surface Stress methodology (CSS) (Gueyffier, Li, Nadim, Scardovelli and Zaleski 1999). A Level-Set formulation of the CSF approach was introduced by Chang et al. (1996). For the second approach, the Ghost Fluid methodology (Liu, Fedkiw and Kang 2000, Kang, Fedkiw and Liu 2000) is most popular, as the alternative Immersed Interface method (Li and Lai 2001) has not yet seen successful im-plementation for flows with discontinuously varying density and viscosity. An advantage of the CSF/CSS methodology is its straightforward implementation, but the disadvantage of regularization is the requirement of very fine grids to accurately approximate the discontinuities. On the other hand, the Ghost-Fluid method allows for very accurate representation of stresses at the interface but its implementation is relatively involved.
A second significant challenge in the implementation of the front-capturing methodology is the representation and advection of the interface. Specifically the accurate representation of interface curvature for the evaluation of surface tension and the volume-conserving advection of fluid phases for substantial in-terface deformation and extended simulation times have been the subject of continuous research in this field (Boersma, Coyajee and Portela 2006).
In the current work, the governing equations for fluid and interface motion are solved numerically on a fixed, uniform Cartesian grid in a three-dimensional rectangular computational domain. The equations for the fluid velocity and in-terface position are integrated sequentially in time: the fluid velocity is advanced given the interface position, and the interface is advanced given the updated flow field. Two different combined Level-Set and Volume-of-Fluid methods are considered for the representation and advection of the interface: The CLSVOF method by Sussman and Puckett (2000) and the MCLS method by Van der Pijl, Segal, Vuik and Wesseling (2005). In these methods, the Level-Set description is particularly convenient for the implementation of jump conditions with the Ghost-Fluid method, while the Volume-of-Fluid methodology is used to ensure the mass-conserving transport of fluid phases.
front-capturing method. Here, our front-front-capturing approach is presented following the traditional approach. In the next chapter, a method is proposed to avoid this numerical coalescence by introducing multiple marker functions for the rep-resentation of separate droplets.
In the following, we first discuss the time integration procedure for the equa-tions of fluid motion using a pressure correction method. For the spatial dis-cretization of the equations of fluid motion, which is discussed in section 3.3, a combined approach is used to deal with discontinuities of the fluid properties at the interface: the viscosity is regularized following the CSF methodology, but surface tension and the discontinuous density are incorporated in the discreti-zation of the Navier-Stokes equations using the Ghost Fluid method. For the selection of the interface representation and advection procedure, two combined Level-Set/Volume-of-Fluid methods are introduced in section 3.4. A discussion of the interface curvature discretization is presented in section 3.5. Subsequently, parallel implementation and performance of the algorithm is discussed in sec-tion 3.6. The chapter concludes with results for several classical problems to evaluate the performance of the single marker front-capturing method.
3.2
Temporal Discretization of the Navier-Stokes
Equations
The Navier-Stokes equations are integrated in time using a pressure correction method (e.g. Van Kan 1986). Equation (2.9) is split into the following predictor and corrector step:
u∗ − un ∆t = − 3 2A (u n) +1 2A u n−1 + 1 ρ 1 ReDimp(u ∗ ) + 1 ReDexp(u n) − Gpn−1 2 + 1 F rz, (3.1) and un+1− u∗ ∆t = − Gp∗ ρ , (3.2) where p∗ = pn+1 2 − pn− 1 2. (3.3)
In (3.1), G represents the discrete gradient and A the discrete advective opera-tor. The advective term is integrated explicitly using the second order Adams-Bashforth method. The diffusive term is split into two parts:
∇ · µ(∇u + ∇uT) = ∇ · (µ∇u) + ∇ · µ∇uT , (3.4)
20 Chapter 3 is defined by: (µu∗ x)x+ µu ∗ y y+ (µu ∗ z)z, (µv∗ x)x+ µv ∗ y y+ (µv ∗ z)z, (µw∗x)x+ µw ∗ y y+ (µw ∗ z)z. (3.5)
Using the incompressibility condition (2.1), the explicit diffusive operator can be written as: − µvyn x− (µw n z)x+ (µvxn)y+ (µwnx)z, µuny x− (µu n x)y− (µwzn)y+ µwny z, (µun z)x+ (µvnz)y− (µunx)z− µvyn z. (3.6)
The discrete operators corresponding to these implicitly and explicitly inte-grated diffusive parts are respectively denoted Dimpand Dexp.
Due to the splitting procedure of the diffusive terms, (3.1) represents three separate linear systems for u∗
, v∗
and w∗
. These systems are solved sequentially with an incomplete Choleski preconditioned conjugate gradient (ICCG) solver (e.g. Golub and Van Loan 1989). By taking the divergence of (3.2) and imposing the incompressibility condition for the velocity at the new time level:
DIV un+1 = 0, (3.7)
a discrete variable coefficient Poisson equation for p∗is obtained:
DIV 1 ρGp ∗ = 1 ∆tDIV (u ∗ ) , (3.8)
where DIV denotes the discrete divergence operator. Equation (3.8) is solved for p∗ with the ICCG method, after which u∗ is corrected toward un+1 using
(3.2). Finally, the pressure is updated according to (3.3).
In (3.1), the material parameters are always derived from the interface po-sition at tn+1 2: µ = µ(φn+1 2), ρ = ρ(φn+ 1 2), (3.9)
where µ(φ) and ρ(φ) are defined by (2.5) and (2.6).
Note that in (3.1), the pressure gradient is evaluated at tn−1
2 and incorpora-tion of jump condiincorpora-tions requires the interface representaincorpora-tion at tn−1
2. To include the correct jump conditions in the discretization of the Poisson problem for p∗
(3.8), the interface representation at both tn+1
2 and tn− 1
2 is required. This is due to the definition of p∗ as the difference between pn+1
2 and pn− 1
2 in (3.3).
Time Step
The time step is determined from an adaptive criterion based on restrictions due to convection, surface tension and gravity:
Here, the time step criterion based on convection is defined by: ∆tc=
h
|u|max+ |v|max+ |w|max
, (3.11)
where h denotes the uniform mesh width. The time step based on surface tension is defined by (Brackbill et al. 1992):
∆ts=
r
W e max(ζ, 1)
4π h
3/2. (3.12)
Following Kang et al. (2000) and Sussman and Puckett (2000), a general time step criterion which includes the effect of gravity is defined by:
∆tf= 2h |un| +p|un|2+ 4hf, (3.13) where f = − 3 2A (u n) +1 2A u n−1 −1ρ 1 Re(D (u n ) − Gpn−12 + 1 F rz . Here D = Dexp+ Dimp denotes the complete discrete diffusive operator. The
latter criterion is based on the consideration of (2.2) by a simplified equation: ut = f, with temporal discretization: un+1 = un+ ∆tf. Requiring a CFL-type
condition for the time step, un+1∆t ≤ h, then leads to the time step criterion
(3.13).
3.3
Spatial Discretization of the Navier-Stokes
Equations
The equations are spatially discretized on a uniform Cartesian grid using finite differences. A standard staggered arrangement of variables is used (Harlow and Welch 1965): Velocity components are located at cell faces, whereas pressure, density, viscosity and Level-Set function are defined at cell centers (figure 3.1). Following Sussman et al. (1994) and Chang et al. (1996), the viscosity is regularized by replacing the Heaviside function H(φ) in (2.5), (2.6) with a con-tinuous representation, denoted Hα(φ):
Hα(φ) = 0 if φ < −α, 1 2(1 + φ α+ 1 πsin( πφ α)) if |φ| ≤ α, 1 if φ > α, (3.14)
where we choose α = 3/2 h, h being the uniform mesh width.
Kang et al. (2000) show that when the viscosity is regularized, the gradients of the velocity are also continuous at the interface and condition (2.10) reduces to:
[p] = κ
22 Chapter 3 z x y v u i,j+1/2,k x ∆ i,j,k+1/2 w y ∆ ∆z pi,j,k φi,j,k ρi,j,k µi,j,k i+1/2,j,k
Figure 3.1: Lay-out of variables in computational cell (i, j, k).
As a result, no jump conditions for the velocity remain and straightforward central second order finite differences are used to approximate derivatives of the velocity, anywhere in the computational domain. For the discretization of the convective and diffusive terms in the Navier-Stokes equations with a central scheme, values of the velocity and viscosity are required at different locations than those available in the staggered arrangement of figure 3.1. Whenever this occurs an average of adjacent values is used, following Harlow and Welch (1965). The interface condition (3.15) represents an explicit jump condition for the pressure which is implemented in the discretization of the pressure gradient with the Ghost-Fluid method (Liu et al. 2000). In the Ghost-Fluid method, known jump conditions are incorporated into a regular finite difference stencil, using the Level-Set function to identify the location of the interface. Given the example in figure 3.2, the pressure gradient at the location of cell face (i +1
2) is approximated as follows: ∂p ∂x i+1 2 =pi+1− pi− [p] ∆x . (3.16)
Taking (3.15) into account, a general formulation for the pressure gradient at cell face (i +1 2) is defined from: ∂p ∂x i+1 2 = pi+1− pi− κΓ/W e ∆x if φi≤ 0 and φi+1> 0, pi+1− pi+ κΓ/W e ∆x if φi> 0 and φi+1≤ 0, pi+1− pi ∆x otherwise. (3.17)
Here κΓ denotes the curvature evaluated at the interface between grid points
i and i + 1. The procedure to determine the interface curvature is discussed in section 3.5.
[p] p p p p i−1 i i+1 i+2
xi−1/2 xi+1/2 xi+3/2
Figure 3.2: One-dimensional example of the discontinuous pressure distribution near an interface (dashed line) with pressure jump [p].
from (3.8). For the one-dimensional example in figure 3.2, the discretization of (3.8) at xi is given by: 1 ρi+1 2 p∗ i+1− p ∗ i − [p ∗ ] ∆x − 1 ρi−1 2 p∗ i − p ∗ i−1 ∆x ! /∆x = 1 ∆t u∗ i+1 2 − u ∗ i−1 2 ∆x , (3.18)
which can also be written as 1 ρi+1 2 p∗ i+1− p ∗ i ∆x − 1 ρi−1 2 p∗ i − p ∗ i−1 ∆x ! /∆x = 1 ∆t u∗ i+1 2 − u ∗ i−1 2 ∆x + 1 ρi+1 2 [p∗ ] ∆x2. (3.19) Because the jump condition for the pressure is incorporated as a source term in the right hand side of the discrete Poisson problem, it can be solved with a fast iterative method. It is easy to see that the matrix associated with the linear system of equations (3.19) is symmetric, which facilitates the use of a conjugate gradient method.
Alternatively, the jump condition for the pressure can be included in the discretization using the Continuous Surface Force (CSF) methodology (Brackbill et al. 1992). Rewriting (3.17) into:
∂p ∂x i+1 2 = pi+1− pi ∆x − κΓ W e H(φi+1) − H(φi) ∆x , (3.20)
the corresponding CSF formulation is found by introducing the smoothed Heav-iside function (3.14): ∂p ∂x i+1 2 =pi+1− pi ∆x − κΓ W e Hα(φi+1) − Hα(φi) ∆x . (3.21)
24 Chapter 3 cell centers, discretization of (3.8) requires values of 1
ρ at cell faces. Following
Liu et al. (2000), the weighted harmonic average is used to obtain values of 1 ρ
at a cell face in the vicinity of the interface, e.g. at xi+1
2: 1 ρi+1 2 = 1 ρiθ + ρi+1(1 − θ) , (3.22)
where ρi= ρ(φi) is directly determined from (2.5) and θ is given by:
θ = |φi| |φi| + |φi+1|
. (3.23)
3.4
Interface Advection
The LS interface representation described in section 2.1 provides all information required for the implementation of interface conditions, e.g. by means of the Ghost Fluid method: the discrete LS function implicitly defines the position of the interface on the computational grid. However, numerical methods for the advection equation of the LS function (2.8) do not conserve the volume of each fluid phase, which adversely affects the quality of long time simulations of incompressible fluid flow. To improve the volume conservation of the fluid phases, the Level-Set method has been combined with the Volume-Of-Fluid method (VOF).
In the VOF method, fluid phases are represented by the discrete volume fraction or Volume-Of-Fluid function. The discrete volume fraction, denoted by ψ, can be defined in terms of the LS function by:
ψi,j,k(t) = 1
∆x∆y∆z Z
Ωi,j,k
H(φ(x, t))dx, (3.24)
where Ω denotes the volume of a computational cell. Consequently, if all fluid in cell (i, j, k) belongs to fluid phase ‘1’, ψi,j,k= 0. If all fluid belongs to phase
‘2’, ψi,j,k = 1. If the cell contains both fluid phase ‘1’ and ‘2’, 0 < ψi,j,k < 1.
Advection of the fluid phases is governed by the following conservation law for the volume fraction:
dψi,j,k dt + 1 ∆x∆y∆z Z ∂Ωi,j,k H(φ)u · n dx = 0, (3.25)
which can also be written as:
ψt+ ∇ · (uψ) = 0. (3.26)
The simplification of (3.25) toward (3.26) is not trivial. The volume fraction function, ψ, only represents the magnitude of the fluid fraction in a compu-tational cell. It does not give explicit information about the distribution of the fluid phases. However, the accuracy of the VOF advection procedure de-pends predominantly on the numerical approximation of the distribution of fluid phases to estimate fluxes of each fluid across cell boundaries. Although this dis-tribution can also be approximated from volume fractions, the LS interface rep-resentation provides a direct approach to identify the location of fluids in a cell, i.e. it facilitates a geometric interface reconstruction for mass-conserving advec-tion of the VOF funcadvec-tion. Therefore, mutual benefits appear when combining the LS and VOF approaches in a single interface advection method.
The first and most established combined LS and VOF method, based on the preliminary work by Bourlioux (1995), is the Coupled Level-Set / Volume-of-Fluid (CLSVOF) method (Sussman and Puckett 2000). Another, more recent approach to combine the LS and VOF methodologies is the Mass Conserving Level-Set (MCLS) interface advection method (Van der Pijl et al. 2005). In either method, both the LS interface and VOF phase representations are ad-vected using conservation equations for the VOF function (3.26) and a similar conservation equation for the LS function:
φt+ ∇ · (uφ) = 0, (3.27)
which is obtained from (2.8) using (2.1)1. VOF advection is treated similarly
in the MCLS and CLSVOF methods, but different discretizations are used to solve the LS equation (3.27). The most distinctive difference between the MCLS and CLSVOF methods is the procedure by which the coupling of LS and VOF functions is established.
A detailed discussion and evaluation of the LS/VOF coupling procedures for the CLSVOF and MCLS methods is presented in chapter 5. In the following, a brief introduction of the CLSVOF and MCLS method is given. The performance of both methods is evaluated by application to a number of advection tests, including the demanding single reversed vortex test.
3.4.1
CLSVOF Method
In this section, the CLSVOF interface advection procedure is described. This procedure consists of two steps. First, the LS and VOF functions are advected individually over a discrete time interval. Subsequently, the LS function is reini-tialized as the signed distance function to piecewise linear interface segments, which are reconstructed from volume fractions, i.e. the advected VOF function values, in interface cells.
1. Advection of the LS and VOF functions: Given un, ψn−1
2 and φn− 1
2,
equations (3.26) and (3.27) are solved using a second order accurate con-servative operator split advection scheme to obtain ψn+1
2 and φn+ 1 2. The
26 Chapter 3 three-dimensional operator split advection scheme is an extension of the two-dimensional implicit-explicit scheme proposed by Puckett, Almgren, Bell, Marcus and Rider (1997). For a general scalar s, it is defined by:
˜ si,j,k = sn− 1 2 i,j,k − F n−1 2 i+1 2,j,k+ F n−1 2 i−1 2,j,k 1 − Un i+1 2,j,k + Un i−1 2,j,k , ˆ si,j,k = ˜ si,j,k− ˜Fi,j+1 2,k+ ˜Fi,j− 1 2,k 1 − Vn i,j+1 2,k+ V n i,j−1 2,k , ¯ si,j,k = ˆ si,j,k− ˆFi,j,k+1 2 + ˆFi,j,k− 1 2 1 − Wn i,j,k+1 2 + Wn i,j,k−1 2 , sn+ 1 2
i,j,k = s¯i,j,k− ˜si,j,k(U n i+1 2,j,k+ U n i−1 2,j,k)− ˆ si,j,k(Vi,j+n 1 2,k+ V n i,j−1 2,k) − ¯si,j,k(W n i,j,k+1 2 + W n i,j,k−1 2). (3.28) Here, Fn− 1 2 i+1 2,j,k = sn− 1 2 i+1 2,j,k Un i+1
2,j,k denotes the flux of s across cell face (i + 1
2, j, k). Velocity U = (U, V, W )
T denotes the CFL-velocity defined
by: U = (u∆t ∆x, v ∆t ∆y, w ∆t ∆z) T. (3.29)
Each timestep, the order of directions in the operator split scheme is in-terchanged to ensure second order accuracy.
The approximation of the Eulerian flux si+1
2,j,k depends on whether s represents the VOF function ψ or LS function φ. The discrete fluxes of ψ are determined from a geometric reconstruction of the fluid phases in each cell (figure 3.3). Assuming a piecewise linear interface in cells where 0 < ψ < 1, the co-advected LS function is used to obtain the normal vector of the reconstructed interface. The discrete volume fraction flux is then defined by:
Z
ΩD
H(φL) dx, (3.30)
where ΩD is the donating region, which contains the fluid that will flow
through a face during a time-step ∆t and φL represents the linear
re-construction of the interface in the cell. An example of a discrete flux, Fi+1
2,j,k, and its donating region is displayed in figure 3.3. To compute the magnitude of the particular VOF flux in figure 3.3, the volume fraction within the donating region of size: u∆t∆y∆z needs to be considered. For additional details about the evaluation of flux F in case of the LS function we refer to Sussman and Puckett (2000).
volume of each fluid phase, i.e.: X i,j,k ψn+ 1 2 i,j,k = X i,j,k ψn− 1 2 i,j,k, (3.31)
the so-called consistency condition defined by: 0 ≤ ψn+
1 2
i,j,k ≤ 1, (3.32)
may not be satisfied. As discussed by Rider and Kothe (1998), flows with large, nonuniform vorticity which severely stretch and deform the interface may produce small under- and overshoots in the VOF function, typically O(10−4). This also occurs when the volume fractions are advected with
the operator split scheme (3.28). Obviously, under- and overshoots result in unphysical values of the volume fraction in a cell. In the CLSVOF method, these are eliminated by truncating values of ψ > 1 toward ψ = 1 and ψ < 0 toward ψ = 0. Due to the truncation, a limited amount of mass transfer between the phases occurs and (3.31) is not exactly maintained. However, results in the next paragraph show that the mass transfer is very small and converges with spatial and temporal refinement.
2. Reinitialization of the LS function: The advection procedure for the LS function yields an interface representation which does not conserve the volume of each fluid phase. Therefore, the advected LS function φn+1
2 is
reinitialized to represent a distance function to the advected VOF function ψn+1
2, which is a volume conserving representation of the fluid phases except for the truncation of under- and overshoots as described in the previous paragraph.
The LS reinitialization is achieved by first reconstructing a piecewise linear interface in each cell for which 0 < ψn+1
2 < 1. The normal vector to a piecewise linear interface segment is obtained from the advected LS function, φn+1
2. The intercept of the linear segment in the cell is obtained from the value of the advected VOF function, ψn+1
2. Finally, in each cell over a width of at least 4 cells at either side of the interface the LS function is assigned to be the signed closest distance to the reconstructed piecewise linear interface.
The CLSVOF advection procedure is formally second order accurate as all procedures including time integration, approximation of fluxes and reconstruc-tion of both the LS and the VOF funcreconstruc-tion are second order accurate. For addi-tional details concerning the CLSVOF advection algorithm we refer to Sussman and Puckett (2000).
28 Chapter 3 u∆t z y x F ∆z y ∆ x+ 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 PSfrag replacemen ts φ ψ
Figure 3.3: Computational cell with flux Fi+1
2,j,k and its donating region u∆t∆y∆z. The interface is indicated by a shaded surface.
Figure 3.4: Graph of the rela-tion between φ and ψ for ∇φ =
1 3 √ 3,1 3 √ 3,1 3 √ 3.
3.4.2
MCLS Method
In this section, the MCLS interface advection procedure is discussed. In the MCLS method, a direct coupling of the LS and VOF functions is pursued in each interface cell using an explicit relation between the LS and VOF func-tions. Starting from an initial LS function, corresponding volume fractions are computed using this LS/VOF function. Similar to the CLSVOF advection procedure, the LS and VOF functions are advected separately. Subsequently, the LS function is locally adjusted so that the volume fractions derived from the explicit LS/VOF relation match the volume fractions of the advected VOF function.
Thus, the MCLS method aims to achieve volume conserving advection of the LS interface representation by advancing the LS and its corresponding VOF function in time with the following scheme:
1. Computation of the VOF function: The volume fractions in each computational cell are directly computed from the LS function. A lin-earization of the LS function around the center of computational cell (i, j, k) is considered. The linearized LS function, denoted φL, is defined by:
φL(x) = φi,j,k+ (x − xi,j,k) · ∇φ i,j,k. (3.33) After substitution of (3.33) into (3.24), the following approximation for the volume fraction in cell (i, j, k) is obtained:
where the derivatives of the LS function are evaluated by central finite differences. The resulting expression, whose compact definition is written:
ψ = f (φ, ∇φ), (3.35)
is equivalent to the analytical relations presented by Scardovelli and Za-leski (2000) and Lorstadt and Fuchs (2004) which define the volume frac-tion in a cube as a funcfrac-tion of the normal vector to the planar interface segment and the distance between cell-corner and plane. A detailed for-mulation of (3.35) is included in appendix A.
2. Advection of the VOF function: Given un, the VOF function ψn−1
2 = f (φn−1 2, ∇φn− 1 2) is advected toward ψn+ 1
2 with the operator split scheme (3.28). The volume fluxes are computed directly from the LS function, using relation (3.35). In the MCLS method, the LS function is advected separately using an unsplit discretization as described in the following paragraph. Therefore, the LS advection does not provide a suitable LS representation of the interface, ˜φ and ˆφ, to compute the intermediate VOF fluxes ˜F and ˆF . Such is the case in the CLSVOF method, which uses the same operator split scheme to advect both the LS and VOF functions. Instead, in the MCLS method intermediate LS functions ˜φ and ˆφ are obtained by applying corrections to respectively φn−1
2 and ˜φ such that ˜
ψ = f ( ˜φ, ∇ ˜φ) and ˆψ = f ( ˆφ, ∇ ˆφ). This correction procedure is described in the fourth part of this summary.
The MCLS method uses the same operator split scheme as the CLSVOF method to advect the VOF function and hence the same problem with under- and overshoots occurs as described in section 3.4.1. Instead of truncating the volume fractions, the MCLS method eliminates the under-and overshoots by volume redistribution, which improves the overall vol-ume conservation of the VOF advection. For further details concerning the volume redistribution procedure, we refer to Van der Pijl (2005). 3. Advection of the LS function: Although interface advection of the
LS function is not strictly volume conserving, it can be used to obtain an approximately volume conserving solution. Given velocity un, (3.27) is
solved to advance the LS function from time level n −1
2 toward a
prelim-inary update of the LS function at time level n +1
2, denoted by φ
∗. After
advection, the LS function is reinitialized to preserve its distance function property. The LS reinitialization in the MCLS method is based on the procedure by Sussman et al. (1994), however modified to improve volume conservation. For details of the LS reinitialization we refer to Van der Pijl (2005).
30 Chapter 3 et al. (2005) argue in favor of a low order method because the diffusion of the low order discretization assures numerical smoothness of the LS function.
4. LS correction: The preliminary advected LS function, φ∗, is corrected
toward φn+1
2 to provide a solution which satisfies:
|ψn+ 1 2 i,j,k − f(φ n+1 2 i,j,k, ∇φ n+1 2 i,j,k)| ≤ , (3.36)
where is a tolerance. In this work = 1 · 10−4.
The essential idea of the correction procedure is the following: Keeping ∇φ at a constant value, the interface equation (3.35) represents an explicit relation between the interface position and the volume fraction in a cell:
φ = g(ψ, ∇φ). (3.37)
Figure 3.4 shows an example of the graph of such a relation between the LS and VOF function for ∇φ = 1
3 √ 3,1 3 √ 3,1 3 √ 3. From figure 3.4 it will be clear that, while assuming a constant gradient of the LS function, the local value of the LS function in a cell can be adjusted such that the required local value of the VOF function is obtained. We note that the local value of the volume fraction is approximated due to the application of a tolerance > 0 in (3.36). Consequently, volume is not exactly conserved in the LS reconstruction procedure.
After the value of the LS function is adjusted locally according to (3.37), the gradient of the corrected LS function will also have altered so that (3.36) does not necessarily hold. Therefore corrections are repeated iter-atively until convergence, as presented by the following scheme:
• φ0= φ∗
; • δψi,j,k = |ψ
n+1 2
i,j,k − f(φ0i,j,k, ∇φ0i,j,k)|;
• Repeat until δψi,j,k< ∀(i, j, k):
φm+1i,j,k = (
g(ψn+ 1 2
i,j,k, ∇φmi,j,k) if δψi,j,k>
φm i,j,k if δψi,j,k≤ ; δψi,j,k = |ψ n+1 2 i,j,k − f(φ m+1 i,j,k, ∇φ m+1 i,j,k)|; • φn+1 2 = φm+1.
The iterative procedure is generally found to converge within 4 to 10 it-erations.
The local application of the corrections is found to have negative impli-cations on interface curvature estimated from the LS function (Coyajee, Herrmann and Boersma 2004, Van der Pijl 2005). A detailed discussion of this phenomenon and improvements of the LS correction procedure are presented in chapter 5.
We remark that the MCLS method appears to be similar to the CLSVOF method in many respects. For both the VOF advection step and the procedure which enforces the coupling of LS and VOF functions, the geometric interface reconstruction is the same: The VOF function is used to define the intercept of a piecewise linear interface segment in mixed cells (0 < ψi,j,k < 1), while the
LS function is used to compute its normal vector. However, formal accuracy of the MCLS method is only first order, against second order for CLSVOF. Although the VOF advection in the MCLS method is potentially second order accurate, discretization of the fluxes in the LS advection equation and the LS correction procedure are only first order accurate, yielding first order accuracy for the entire method.
3.4.3
Results of 2D Advection Test Cases
In this section, results of the MCLS and CLSVOF interface advection meth-ods are presented, applied to several two-dimensional test cases. In the past, Van der Pijl (2005) and Sussman and Puckett (2000) reported results for Zale-sak’s rotation of a notched disk, however using different error definitions. Here, the MCLS and CLSVOF methods are directly compared for Zalesak’s test case. In addition, first results of the MCLS and CLSVOF methods are presented for the reversed single vortex problem introduced by Rider and Kothe (1998). Zalesak’s Rotating Slotted Disk Test
Zalesak’s test case for the advection of a notched disk (Zalesak 1979) has be-come a standard test case in literature. A segment of fluid in the shape of a notched disk is placed in a rectangular computational domain of dimensions L × L (figure 3.5). The center of the disk of radius R0 = 203L is located at
position (x0, y0) and the dimension of the width of notch and upper bridge is
w = 1
3R0. The fluid segment is passively advected, either with the CLSVOF
or the MCLS method, using a velocity field representing a rotational flow field around the center of the computational domain:
u = (1
2L − y), v = (x − 1
2L). (3.38)
The advection time is chosen such that the fluid segment completes exactly one revolution in the analytical solution.
To compare results of the CLSVOF method with those of the MCLS method (Van der Pijl et al. 2005), we use (x0, y0) = (12L,34L) and ∆t = 6.28Lπh . Also in
32 Chapter 3 0 L w w R 0 y x0 L x y
Figure 3.5: Initial configuration for Zalesak’s advection test.
interface length, l(φ), defined by: l(φ) =
Z
Ω
δ(φ)|∇φ| dΩ, (3.39)
and the shift of the interface by means of error norm ||e||2, defined:
||e||2= R Ω φ−φ0 L 2 δ(φ)|∇φ| dΩ R Ωδ(φ)|∇φ| dΩ 1 2 , (3.40)
with φ0 the exact, analytical solution for the Level-Set function. In (3.39)
and (3.40), integrals are evaluated using the midpoint rule, the gradients are approximated by central second order finite differences and the delta function is replaced by a regularized version, denoted δα:
δα(x) = 1 2α 1 + cosπx α if |x| ≤ α, 0 otherwise, (3.41)
where α = 3/2h, h being the uniform mesh width. In figure 3.6, results are presented for the relative error in interface length, |l(φ)/l(φ0
10-3 10-2 10-1 10-3 10-2 h || e|| 2 10-3 10-2 10-1 10-2 10-1 100 h |l (φ )/l (φ 0) − 1|
Figure 3.6: Errors for interface shift (left) and interface length (right) for Za-lesak’s advection test. Comparison between results of MCLS (squares, dashed line) and CLSVOF (circles, solid line) advection methods.
appears to be just slightly more accurate than the MCLS method, for both error measurements.2
Interface contours of the CLSVOF method are shown in figure 3.7 for com-parison with figure 8 in Van der Pijl et al. (2005). At the highest resolution no notable difference between contours can be observed. At the lowest resolution, the CLSVOF result is clearly more accurate, not showing the change in topology as observed for the MCLS method.
While the CLSVOF is formally second order and the MCLS method first order accurate, the difference in accuracy between the CLSVOF and MCLS methods for Zalesak’s test case is small. On the one hand the accuracy of the CLSVOF results is limited due to the nature of Zalesak’s test case: The sharp corners of the slotted disk can only be resolved with first order accuracy. On the other hand the MCLS method is expected to perform well resolving the sharp corners due to the local LS reconstruction procedure.
Fluid Disk in a Single Reversed Vortex
In this section, the MCLS and CLSVOF advection methods are compared for the deformation of a fluid disk in a single reversed vortex, a test introduced by
2In (Van der Pijl et al. 2005), the exact interface length, l(φ0), is defined incorrectly. It
should read: l(φ0 ) = 2R0−w+ 2 s R02− „ 1 2w «2 + 2πR0−2 arcsin „ w 2R0 « R0. (3.42)
34 Chapter 3 50 x 50 0.3 0.4 0.5 0.6 0.7 0.6 0.7 0.8 0.9 100 x 100 0.3 0.4 0.5 0.6 0.7 0.6 0.7 0.8 0.9 150 x 150 0.3 0.4 0.5 0.6 0.7 0.6 0.7 0.8 0.9 200 x 200 0.3 0.4 0.5 0.6 0.7 0.6 0.7 0.8 0.9
Figure 3.7: Comparison of analytical and numerical (CLSVOF) solutions for Zalesak’s advection test. Results displayed for grids of 502, 1002, 1502and 2002
computational cells.
Rider and Kothe (1998). The results are also compared to those presented by Scardovelli and Zaleski (2003).
A fluid disk of radius 3
20 is placed at point ( 1 2,
3
4) in a unit square domain.
The solenoidal velocity field u = (u, v) is prescribed by: u = − sin2(πx) sin(2πy) cos πt T
, v = sin(2πx) sin2(πy) cos πt
T
.
(3.43)
0 0.25 0.5 0.75 1 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1 0.25 0.5 0.75 1
Figure 3.8: Fluid disk in a single reversed vortex: Interface configuration and velocity field at t = 0 (left) and t = T /2 (right). Result at t = T /2 computed with the CLSVOF method on a grid of 1282cells.
returns to its initial position. The accuracy of the MCLS and CLSVOF methods is assessed by comparing the initial and final position of the fluid after a single period T = 2 on grids containing 322, 642, 1282 and 2562 computational cells.
The error of the discrete volume fractions between the exact and computed solutions is defined by:
Ea= h2
X
i,j
|ψi,j(t = T ) − ψi,j(t = 0)|, (3.44)
where h is the uniform mesh width and ψi,j denotes the volume fraction in cell
(i, j). The error of the relative mass/volume which is lost or gained in respect to the original mass/volume is defined by:
Em= X i,j ψi,j(t = T ) − ψi,j(t = 0) ψi,j(t = 0) . (3.45)
In table 3.1, we show that the CLSVOF method is more accurate than the MCLS method for the single reversed vortex test. In figure 3.9, a comparison between initial and final interface contours is made both for the CLSVOF and the MCLS method on 322 grid. Not only a difference in the magnitude of
