Stéphane Popinet
email: s.popinet@niwa.co.nz

An accurate adaptive solver for surface-tension-driven interfacial flows

Keywords: 65M50, 76D05, 76D25, 76D45, 76M12, 76T10 Adaptive mesh refinement, Surface tension, Navier-Stokes equations, Height-Function, Parasitic currents, Volume-Of-Fluid, Octree Numerical scheme

A method combining an adaptive quad/octree spatial discretisation, geometrical Volume-Of-Fluid interface representation, balanced-force continuum-surface-force surface tension formulation and height-function curvature estimation is presented. The extension of these methods to the quad/octree discretisation allows adaptive variable resolution along the interface and is described in detail. The method is shown to recover exact equilibrium (to machine accuracy) between surfacetension and pressure gradient in the case of a stationary droplet, irrespective of viscosity and spatial resolution. Accurate solutions are obtained for the classical test case of capillary wave oscillations. An application to the capillary breakup of a jet of water in air further illustrates the accuracy and efficiency of the method. The source code of the implementation is freely available as part of the Gerris flow solver.

Introduction

Flows of immiscible fluids are ubiquitous in Nature: waves on the sea, waterfalls, rain drops or a garden sprinkler are familiar examples. In general air-water flows exhibit the typical features associated with many important two-phase-flow phenomena: relatively large density ratios, high surface-tension and low viscosity at practical length scales. Furthermore most of these flows tend to generate complex and evolving interface geometries on spatial scales ranging over several orders of magnitude. This is true not only of obviously complex phenomena such as wave breaking or jet atomisation, but also of the apparently simpler sub-phenomena such as jet/droplet breakup and coalescence.

Numerical methods using an implicit representation of the interface such as Volume-Of-Fluid (vof) [START_REF] Hirt | Volume of fluid method for the dynamics of free boundaries[END_REF][START_REF] Li | Calcul d'interface affine par morceaux[END_REF][START_REF] Rider | Reconstructing volume tracking[END_REF][START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF] or Levelset (ls) [START_REF] Sussman | A level set approach for computing solutions to incompressible twophase flow[END_REF][START_REF] Sussman | An adaptive level set approach for incompressible two-phase flows[END_REF][START_REF] Sochnikov | Level set calculations of the evolution of boundaries on a dynamically adaptive grid[END_REF][START_REF] Kohno | Numerical analysis of moving interfaces a level set method coupled with adaptive mesh refinement[END_REF][START_REF] Zheng | Adaptive unstructured remeshing-II: Application to two-and three-dimensional level-set simulations of multiphase flow[END_REF][START_REF] Liu | Sharp interface cartesian grid method ii: A technique for simulating droplet interactions with surfaces of arbitrary shape[END_REF][START_REF] Marchandise | A stabilized finite element method using a discontinuous level set approach for solving two phase incompressible flows[END_REF][START_REF] Marchandise | A stabilized finite element method using a discontinuous level set approach for the computation of bubble dynamics[END_REF][START_REF] Desjardins | An accurate conservative level set/ghost fluid method simulating turbulent atomization[END_REF] can robustly and efficiently represent evolving, topologically complex interfaces. The accurate representation of surface tension within these methods is typically more delicate than for schemes using an explicit representation of the interface [START_REF] Fyfe | Surface tension and viscosity with Lagrangian hydrodynamics on a triangular mesh[END_REF][START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF][START_REF] Shin | Accurate representation of surface tension using the level contour reconstruction method[END_REF][START_REF] Dai | Adaptive tetrahedral meshing in free-surface flow[END_REF] and care must be taken to avoid any imbalance between discrete surface tension and pressure gradient terms [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF][START_REF] Gueyffier | Volume of fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF][START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF][START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF][START_REF] Renardy | PROST -A parabolic reconstruction of surface tension for the volume-offluid method[END_REF][START_REF] Harvie | An analysis of parasitic current generation in volume of fluid simulations[END_REF][START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF]. For vof methods the accurate evaluation of the geometrical properties of the interface such as curvature is also an active research topic [START_REF] Ginzburg | Two-phase flows on interface refined grids modeled with VOF, staggered finite volumes, and spline interpolants[END_REF][START_REF] Renardy | PROST -A parabolic reconstruction of surface tension for the volume-offluid method[END_REF][START_REF] Cummins | Estimating curvature from volume fractions[END_REF][START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF][START_REF] Sussman | Improvements for calculating two-phase bubble and drop motion using an adaptive sharp interface method[END_REF]. Finally dynamic mesh-adaptive methods can deal efficiently with phenomena involving a wide range of spatial scales [START_REF] Almgren | A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations[END_REF][START_REF] Jeong | Finite element analysis of transient fluid flow with free surface using VOF (volumeof-fluid) method and adaptive grid[END_REF][START_REF] Sussman | A parallelized, adaptive algorithm for multiphase flows in general geometries[END_REF][START_REF] Sussman | An adaptive level set approach for incompressible two-phase flows[END_REF][START_REF] Sochnikov | Level set calculations of the evolution of boundaries on a dynamically adaptive grid[END_REF][START_REF] Greaves | A quadtree adaptive method for simulating fluid flows with moving interfaces[END_REF][START_REF] Kohno | Numerical analysis of moving interfaces a level set method coupled with adaptive mesh refinement[END_REF][START_REF] Theodorakakos | Simulation of sharp gas-liquid interface using VOF method and adaptive grid local refinement around the interface[END_REF][START_REF] Dai | Adaptive tetrahedral meshing in free-surface flow[END_REF][START_REF] Zheng | Adaptive unstructured remeshing-II: Application to two-and three-dimensional level-set simulations of multiphase flow[END_REF][START_REF] Losasso | Spatially adaptive techniques for level set methods and incompressible flow[END_REF][START_REF] Min | A second order accurate level set method on non-graded adaptive Cartesian grids[END_REF][START_REF] Compère | Transient adaptivity applied to two-phase incompressible flows[END_REF].

In the following I will describe a new combination of a vof and quad/octree adaptive mesh refinement method. The basis for the scheme is the adaptive, incompressible quad/octree Gerris solver described in [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. It is coupled with a classical geometrical vof scheme which is generalised to work on the quad/octree in order to allow variable spatial resolution along the interface. This is in contrast with many previous implementations of adaptive methods with vof or ls which restricted adaptivity to regions away from the interface (i.e. resolution along the interface was kept constant) [START_REF] Jeong | Finite element analysis of transient fluid flow with free surface using VOF (volumeof-fluid) method and adaptive grid[END_REF][START_REF] Sochnikov | Level set calculations of the evolution of boundaries on a dynamically adaptive grid[END_REF][START_REF] Wang | Finite-volume-type VOF method on dynamically adaptive quadtree grids[END_REF][START_REF] Greaves | A quadtree adaptive method for simulating fluid flows with moving interfaces[END_REF][START_REF] Kohno | Numerical analysis of moving interfaces a level set method coupled with adaptive mesh refinement[END_REF][START_REF] Theodorakakos | Simulation of sharp gas-liquid interface using VOF method and adaptive grid local refinement around the interface[END_REF][START_REF] Losasso | Simulating water and smoke with an octree data structure[END_REF][START_REF] Zheng | Adaptive unstructured remeshing-II: Application to two-and three-dimensional level-set simulations of multiphase flow[END_REF][START_REF] Greaves | Simulation of viscous water column collapse using adapting hierarchical grids[END_REF][START_REF] Nikolopoulos | Three-dimensional numerical investigation of a droplet impinging normally onto a wall film[END_REF][START_REF] Min | A second order accurate level set method on non-graded adaptive Cartesian grids[END_REF][START_REF] Herrmann | A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids[END_REF]. A recent article by Malik et al [START_REF] Malik | Adaptive VOF with curvature-based refinement[END_REF] describes a two-dimensional quadtree implementation allowing variable resolution along the interface but their study is limited to the advection problem (not coupled to the momentum equation). Block-structured amr and levelset methods were used in [START_REF] Sussman | An adaptive level set approach for incompressible two-phase flows[END_REF][START_REF] Sussman | A parallelized, adaptive algorithm for multiphase flows in general geometries[END_REF] to provide variable resolution along the interface. Here we present a novel quad/octree discretisation.

Surface tension representation draws on the "balanced-forced" continuum surface force (csf) concept introduced by Renardy and Renardy [START_REF] Renardy | PROST -A parabolic reconstruction of surface tension for the volume-offluid method[END_REF] and Francois et al [START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF], coupled with a height-function (hf) curvature calculation [START_REF] Rudman | Volume-tracking methods for interfacial flows calculations[END_REF][START_REF] Cummins | Estimating curvature from volume fractions[END_REF] generalised for the quad/octree. The classical hf curvature calculation becomes inconsistent at low interface resolution. This is addressed by a generalised version of the hf curvature concept.

The robustness and accuracy of the generalised surface-tension method is assessed for classical test cases: curvature estimates, stationary droplet with surface tension, capillary waves and rising bubbles and compared with published results from other methods. Particular attention is paid to the issue of "parasitic currents" around a stationary droplet [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF][START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF][START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF][START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF][START_REF] Sussman | Improvements for calculating two-phase bubble and drop motion using an adaptive sharp interface method[END_REF][START_REF] Herrmann | A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids[END_REF], and we improve upon existing methods by enabling a unique combination of varying resolution along the interface, improved height-function curvature discretisation and a balanced-force approach. Finally the method is applied to the adaptive solution of the surface-tension-driven breakup of a three-dimensional cylindrical jet of water in air.

The full implementation of the method as well as most of the test cases presented in this article are available under the General Public License on the Gerris web site [START_REF] Popinet | The Gerris Flow Solver[END_REF][START_REF] Popinet | Gerris Test Suite[END_REF].

Numerical scheme

The numerical scheme is a direct extension of the scheme described in [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. Consequently I will only give a summary of the elements of the method which are not specific to the Volume-Of-Fluid and surface tension implementations. Details regarding the quad/octree spatial discretisation, multilevel Poisson solver, approximate projection and momentum advection scheme can be found in [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF].

Temporal discretisation

The incompressible, variable-density, Navier-Stokes equations with surface tension can be written

ρ (∂ t u + u • ∇u) = -∇p + ∇ • (2 µ D) + σ κ δ s n, ∂ t ρ + ∇ • (ρ u) = 0, ∇ • u = 0,
with u = (u, v, w) the fluid velocity, ρ ≡ ρ(x, t) the fluid density, µ ≡ µ(x, t) the dynamic viscosity and D the deformation tensor defined as D ij ≡ (∂ i u j + ∂ j u i )/2. The Dirac distribution function δ s expresses the fact that the surface tension term is concentrated on the interface; σ is the surface tension coefficient, κ and n the curvature and normal to the interface.

For two-phase flows we introduce the volume fraction c(x, t) of the first fluid and define the density and viscosity as

ρ(c ˜) ≡ c ˜ρ1 + (1 -c ˜) ρ 2 , (1) µ
(c ˜) ≡ c ˜µ1 + (1 -c ˜) µ 2, (2)
with ρ 1 , ρ 2 and µ 1 , µ 2 the densities and viscosities of the first and second fluids respectively. Field c ˜is either identical to c or is constructed by applying a smoothing spatial filter to c. Using a smoothed field to define the viscosity was found to improve the results for some of the test cases I studied. Unless otherwise specified spatial filtering is not used in the applications presented in the following. When spatial filtering is used, field c ˜is constructed by averaging the four (respectively eight in 3d) corner values of c obtained by bilinear interpolation from the cell-centred values. When spatial filtering is used, the properties associated with the interface are thus "smeared" over three discretisation cells.

The advection equation for the density can then be replaced with an equivalent advection equation for the volume fraction

∂ t c + ∇ • (c u) = 0.
A staggered in time discretisation of the volume-fraction/density and pressure leads to the following formally second-order accurate time discretisation

ρ n+ 1 2 u n+1 -u n ∆t + u n+ 1 2 • ∇u n+ 1 2 = -∇p n+ 1 2 + ∇ • µ n+ 1 2 (D n + D n+1 ) + (σ κ δ s n) n+ 1 2 , c n+ 1 2 -c n-1 2 ∆t + ∇ • (c n u n ) = 0, ∇ • u n = 0
This system is further simplified using a classical time-splitting projection method [START_REF] Chorin | On the convergence of discrete approximations to the Navier-Stokes equations[END_REF] ρ

n+ 1 2 u ⋆ -u n ∆t + u n+ 1 2 • ∇u n+ 1 2 = ∇ • µ n+ 1 2 (D n + D ⋆ ) + (σ κ δ s n) n+ 1 2 , (3)
c n+ 1 2 -c n-1 2 ∆t + ∇ • (c n u n ) = 0, (4)
u n+1 = u ⋆ - ∆t ρ n+ 1 2 ∇p n+ 1 2 , (5)
∇ • u n+1 = 0
which requires the solution of the Poisson equation

∇ • ∆t ρ n+ 1 2 ∇p n+ 1 2 = ∇ • u ⋆ ( 6
)
This equation is solved using the quad/octree-based multilevel solver described in detail in [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. The iterative solution procedure is stopped whenever the maximum relative change in the volume of any discretisation element (due to the remaining divergence of the velocity field) is less than a given threshold γ p . Unless otherwise specified this threshold is set to 10 -3 . The number of Gauss-Seidel relaxations per level is set to four.

It is well-known that the standard multigrid scheme can exhibit slow convergence in the case of elliptic equations with discontinuous coefficients and/or source terms [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF][START_REF] Tatebe | The multigrid preconditioned conjugate gradient method[END_REF][START_REF] Chan | Robust multigrid methods for nonsmooth coefficient elliptic linear systems[END_REF].

Depending on the problem and interface topology, this can lead to an important degradation of performance for large-density-ratio flows. While this was not an issue for the applications presented in this paper, it is an aspect of the method which warrants further improvements.

The discretised momentum equation ( 3) can be re-organised as

ρ n+ 1 2 ∆t u ⋆ -∇ • µ n+ 1 2 D ⋆ = ∇ • µ n+ 1 2 D n + (σ κ δ s n) n+ 1 2 + ρ n+ 1 2 u n ∆t - u n+ 1 2 • ∇u n+ 1 2 (7)
where the right-hand-side depends only on values at time n and n + 1/2. This is an Helmholtz-type equation which can be solved using a variant of the multilevel Poisson solver.

The resulting Crank-Nicholson discretisation of the viscous terms is formally second-order accurate and unconditionally stable. The criterion for convergence of the multilevel solver is the relative error γ u in each component of the velocity field. Unless otherwise specified this threshold is set to 10 -6 . The number of Gauss-Seidel relaxations per level is set to four.

The velocity advection term u n+ 1 2

• ∇u n+ 1 2 is estimated using the Bell-Colella-Glaz secondorder unsplit upwind scheme [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF][START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. This scheme is stable for cfl numbers smaller than one.

Spatial discretisation

Space is discretised using a graded quadtree partitioning (octree in three dimensions). We refer the reader to [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] and references therein for a more detailed presentation of this data structure and just give a summary of the definitions necessary for the description of the algorithms presented in this article:

Root cell. The base of the cell tree. The root cell does not have a parent cell and all cells in the tree are descendants of the root cell.

Children cells. The direct descendants of a cell. Cells other than leaf cells have four children in two dimensions (quadtree) and eight in three dimensions (octree).

Parent cell.

The direct ancestor of a given cell.

Leaf cells.

The highest cells in the cell tree. Leaf cells do not have children.

Cell level.

By convention the root cell is at level zero and each successive generation increases the cell level by one.

Coarser cell.

Cell A is coarser than cell B if level(A) < level(B) and conversely for finer cells.

All the variables (components of the momentum, pressure and passive tracers) are collocated at the centre of each square in 2d (resp. cubic in 3d) discretisation volume. Consistently with a finite-volume formulation, the variables are interpreted as the volume-averaged values for the corresponding discretisation volume. The choice of a collocated definition of all variables makes momentum conservation simpler when dealing with mesh adaptation [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF].

It is also a necessary choice in order to use the Godunov momentum advection scheme of Bell, Colella and Glaz [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF], and it simplifies the implementation of the Crank-Nicholson discretisation of the viscous terms; however one has to be careful to avoid the classic problem of decoupling of the pressure and velocity field.

To do so, an approximate projection method [START_REF] Almgren | Approximate projection methods: part I. inviscid analysis[END_REF][START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] is used for the spatial discretisation of the pressure correction equation ( 5) and the associated divergence in the Poisson equation [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF]. In a first step the auxiliary cell-centred velocity field u ⋆ c is computed using equation ( 7). An auxiliary face-centred velocity field u ⋆ f is then computed using averaging of the cellcentred values on all the faces of the Cartesian discretisation volumes. When faces are at the boundary between different levels of refinement of the quad/octree mesh, averaging is performed so as to guarantee consistency of the corresponding volume fluxes (see [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] for details).

The divergence of the auxiliary velocity field appearing on the right-hand-side of equation ( 6) is then computed for each control volume as the finite-volume approximation

∇ • u ⋆ = 1 ∆ f u ⋆ f • n f ,
with n f the unit normal vector to the face and ∆ the length scale of the control volume.

After solving equation ( 6), the pressure correction is applied to the face-centred auxiliary field

u n+1 f = u ⋆ f - ∆t ρ n+ 1 2 f ∇ f p n+ 1 2 , (8)
where ∇ f is a simple face-centred gradient operator (consistent at coarse/fine volume boundaries, see section 4.1 of [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]). The resulting face-centred velocity field u n+1 f is exactly non-divergent by construction.

The cell-centred velocity field at time n + 1 is obtained by applying a cell-centred pressure correction

u n+1 c = u ⋆ c - ∆t ρ n+ 1 2 f ∇ f p n+ 1 2 c , (9)
where the || c operator denotes averaging over all the faces delimiting the control volume. The resulting cell-centred velocity field u n+1 c is approximately divergence-free.

The Poisson equation [START_REF] Sussman | An adaptive level set approach for incompressible two-phase flows[END_REF], advection-diffusion equation [START_REF] Sochnikov | Level set calculations of the evolution of boundaries on a dynamically adaptive grid[END_REF] and pressure corrections [START_REF] Kohno | Numerical analysis of moving interfaces a level set method coupled with adaptive mesh refinement[END_REF] and ( 9) require an estimate of the face-centred values of the density ρ n+1/2 f or viscosity µ n+1/2 f . Both quantities are obtained using the volume-fraction-weighted averages (1) and (2) which in turn requires an estimate of c n+1/2 f . This face-centred volume fraction is estimated using a simple average of the cell-centred values in the case of neighbouring cells at the same level. At fine-coarse and coarse-fine boundaries a second-order interpolation stencil is used, identical to the stencil used to evaluate the face-centred gradient operator ∇ f (see section 4.1 of [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]). Other choices are possible, including a local face-centred vof reconstruction which would eliminate the smoothing effect of the interpolation stencil at coarse-fine boundaries. Preliminary tests suggest that the interpolation procedure works well in practice however.

Volume Of Fluid advection scheme

To solve the advection equation ( 4) for the volume fraction we use a piecewise-linear geometrical Volume Of Fluid (vof) scheme [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF] generalised for the quad/octree spatial discretisation. Geometrical vof schemes classically proceed in two steps:

1. Interface reconstruction 2. Geometrical flux computation and interface advection

Interface reconstruction

In a piecewise-linear vof scheme the interface is represented in each cell by a line (resp. plane in three dimensions) described by equation

m • x = α, ( 10
)
where m is the local normal to the interface and x is the position vector. Given m and the local volume fraction c, α is uniquely determined by ensuring that the volume of fluid contained in the cell and lying below the plane is equal to c. This volume can be computed relatively easily by taking into account the different ways a square (resp. cubic) cell can be cut by a line (resp. plane) which leads to matched linear and quadratic (resp. cubic) functions of α. This step has been described in detail in several papers [START_REF] Gueyffier | Volume of fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF][START_REF] Scardovelli | Analytical relations connecting linear interfaces and volume fractions in rectangular grids[END_REF]. In what follows we will assume that given a normal m, a volume fraction c and coefficient α of equation [START_REF] Liu | Sharp interface cartesian grid method ii: A technique for simulating droplet interactions with surfaces of arbitrary shape[END_REF] defined in a coordinate system with origin at the cell centre and in which the cell size is unity, we have the relations:

c = V(m, α), α = V -1 (m, c).
In practice, we use the V and V -1 routines implemented by Ruben Scardovelli, based on analytical formulae for the solution of quadratic (resp. cubic) equations [START_REF] Scardovelli | Analytical relations connecting linear interfaces and volume fractions in rectangular grids[END_REF]. The routines have been well tested and are free from inconsistencies (such as V(m, V -1 (m, c)) c) which may occur due to round-off errors in limiting cases.

A number of schemes have been proposed for interface normal estimation [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF][START_REF] Scardovelli | Interface reconstruction with least-square fit and split eulerian-lagrangian advection[END_REF][START_REF] Aulisa | Interface reconstruction with least-squares fit and split advection in three-dimensional Cartesian geometry[END_REF]. Most of these schemes only require information in a compact neighbourhood of the cell considered: typically a 3 × 3 (resp. 3 × 3 × 3) regular Cartesian stencil. Given this stencil, the methods use finite-difference estimates and/or minimisation techniques and return an estimate for m. Generalising this approach to a quad/octree discretisation is trivial if the interfacial region (i.e. the compact neighbourhood) is resolved with a constant resolution. In this case the local discretisation simply reduces to a regular Cartesian discretisation [START_REF] Greaves | A quadtree adaptive method for simulating fluid flows with moving interfaces[END_REF][START_REF] Theodorakakos | Simulation of sharp gas-liquid interface using VOF method and adaptive grid local refinement around the interface[END_REF][START_REF] Wang | Finite-volume-type VOF method on dynamically adaptive quadtree grids[END_REF].

In order to generalise the method to variable spatial resolution along the interface, we need to be able to reconstruct a local regular Cartesian stencil when cells in the compact neighbourhood vary in size. For the sake of clarity we describe the stencil-computation algorithm for two-dimensional quadtrees. The extension to three dimensions is straightforward.

Given a reference cell C centred on x 0 , of size ∆, and containing the interface (0 < c < 1), the stencil-computation algorithm proceeds as follows:

Algorithm 1 Stencil computation (C) 1. set X (1, 0), Y (0, 1).
2. for each position x i,j

x 0 + ∆ (i X + j Y), -1 i 1, -1 j 1 a) locate the smallest cell N of size larger than or equal to ∆ containing x i, j b) set the stencil volume fraction c i, j to the volume fraction of a (virtual) cell N ∆ of size ∆, centred on x i, j and entirely contained in N

Step 1.2.a of the algorithm is performed using a conditional traversal of the quadtree starting from the root cell and stopping whenever either a leaf cell is reached or the size of the cell is equal to ∆. This leads to an efficient point-location algorithm with cost of order O(log n)

where n is the number of leaf cells in the quadtree.

For step 1.2.b, several cases must be considered:

Algorithm 2 Equivalent volume fraction (C,x i, j ,N)
1. the size of N is equal to the size of C, then N ∆ = N and c i,j is simply set to the value of the volume fraction in cell N 2. the size of N is larger than ∆ (i.e. N is a coarse neighbour of C and is also a leaf cell) a) N does not contain the interface: c i, j is set to the volume fraction in cell N ( 0 or 1) b) N contains the interface: by definition N ∆ does not exist in this case and the equivalent volume fraction c i, j needs to be computed from the interface reconstructed in cell N Note that in case 2.1, cell N is not necessarily a leaf cell, however consistent values for the volume fraction can easily be defined on all levels of the quadtree. This is done before interface reconstruction by traversing the quadtree from the leaf cells to the root cell and setting the volume fractions of non-leaf cells to the average of their children's volume fractions.

For the treatment of case 2.2.b, we will assume that the equation m • x = α N is known for the fragment of interface contained in cell N . In order to be able to compute the equivalent volume fraction of the virtual cell N ∆ using function V, we need to transform the interface equation into the coordinate system centred on N ∆ . The coordinate transformation is simply

x ′ = n ∆ x + x N -x i,j ∆ , (11)
where n ∆ and x N are the size and position of cell N . The interface equation in the coordinate system of the virtual cell N ∆ is then

m • x ′ = n α N + m ∆ • (x N -x i,j ), (12)
and the corresponding stencil volume fraction is given by

c i,j = V(m, n α N + m ∆ • (x N -x i,j )). (13)
Note that in practice we restrict the quad/octree discretisation so that the sizes of neighbouring cells never differ by more than a factor of two [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] i.e. n is always equal to two in the expressions above.

Once the c i, j stencil has been filled, a standard Cartesian routine for computing the interface normal m is called and the interface is reconstructed by evaluating α = V -1 (m, c 0,0 ). In practice we use the Mixed-Youngs-Centred (myc) implementation of Aulisa et al [START_REF] Aulisa | Interface reconstruction with least-squares fit and split advection in three-dimensional Cartesian geometry[END_REF] to evaluate m.

As should be clear from the description of algorithm 2, reconstructing an interface in a cell will require access to the reconstructed interface in any coarser neighbouring cell, however access to cells finer than the cell considered is never required. Taking this into account, the general interface reconstruction algorithm can be summarised as:

Algorithm 3 Interface reconstruction 1. for each non-leaf cell (traversing from leaf to root) a) set the volume fraction as the average of the children's volume fractions 2. for each cell C containing the interface (traversing from root to leaf) a) fill the c i,j regular Cartesian stencil using Algorithm 1 b) compute m using c i,j and the myc scheme c) compute α = V -1 (m, c 0,0 ) d ) store m and α as state variables of cell C

Interface advection and geometrical flux computation

The basis for volume fraction advection is the direction-split scheme originally implemented by DeBar [START_REF] Debar | Fundamentals of the KRAKEN code[END_REF] and re-implemented by others [START_REF] Rudman | Volume-tracking methods for interfacial flows calculations[END_REF][START_REF] Rider | Reconstructing volume tracking[END_REF][START_REF] Lörstad | High-order surface tension VOF-model for 3D bubble flows with high density ratio[END_REF]. Considering a single two-dimensional Cartesian cell, one timestep of this scheme can be summarised as

c ⋆ i,j V ⋆ i,j = c n-1/2 i, j V n-1/2 i,j + φ n-1/2 i-1/2,j -φ n-1/2 i+1/2, j , (14)
c n+1/2 i,j V n+1/2 i, j = c ⋆ i, j V ⋆ i,j + φ ⋆ i,j -1/2 -φ ⋆ i,j +1/2 , (15)
where V i, j is the effective volume of the cell and φ i+1/2, j is the volume flux of the first phase through the right-hand boundary of the cell (respectively left-hand, top and bottom boundaries for the other indices). The effective volume is defined by taking into account the compression or expansion of the cell due to the divergent one-dimensional velocity field, which gives

V ⋆ i, j ≡ V n-1/2 i,j + (u n i-1/2, j -u n i+1/2,j ) ∆ ∆t, V n+1/2 i,j ≡ V ⋆ i,j + (u n i, j -1/2 -u n i, j +1/2
) ∆ ∆t.

Since the face-centred velocity field u n f is discretely non-divergent, the second relation

reduces to V n+1/2 i, j ≡ V n-1/2 i,j
≡ ∆ 2 . Also, in practice the order of the split advection directions is reversed at each timestep. In two dimensions, equations ( 14) and ( 15) are equivalent to the Eulerian-implicit and Eulerian-explicit formulations of [START_REF] Scardovelli | Interface reconstruction with least-square fit and split eulerian-lagrangian advection[END_REF] respectively and the scheme as a whole is identical to the Eulerian-implicit-explicit scheme described by [START_REF] Rider | Reconstructing volume tracking[END_REF]. Extension of this scheme to three dimensions is straightforward.

Several schemes have been proposed for the computation of the volume fluxes φ. Classical analytical schemes can be used in principle but usually lead to numerical diffusion of the interface. The original Volume-Of-Fluid scheme of Hirt and Nichols [START_REF] Hirt | Volume of fluid method for the dynamics of free boundaries[END_REF] relied on a combination of diffusive and anti-diffusive advection schemes to limit interface diffusion, an idea initially proposed by Zalesak [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF] and subsequently expanded by others [START_REF] Rudman | Volume-tracking methods for interfacial flows calculations[END_REF][START_REF] Ubbink | A method for capturing sharp fluid interfaces on arbitrary meshes[END_REF]. In contrast to this analytical approach, volume fluxes can also be estimated using the geometry of the reconstructed interface [START_REF] Debar | Fundamentals of the KRAKEN code[END_REF][START_REF] Noh | SLIC/simple line interface calculation[END_REF][START_REF] Li | Calcul d'interface affine par morceaux[END_REF]. This geometrical approach is efficient and simple to implement for Cartesian discretisation elements but is not easily adapted to more general spatial discretisations. As will be shown below this is not an issue when using a quad/octree spatial discretisation. The principle of geometrical flux estimation is illustrated in Figure 1. The total volume which will be fluxed to the right-hand neighbour is delimited with a dashed line. The fraction of this volume occupied by the first phase is indicated by the dark grey triangle. By definition the area of this triangle is an estimate of the volume flux φ i+1/2,j . This area can easily be calculated using an approach similar to that presented in the preceding section (relations [START_REF] Marchandise | A stabilized finite element method using a discontinuous level set approach for solving two phase incompressible flows[END_REF], ( 12) and ( 13)). Generalising this advection scheme to a quad/octree discretisation only involves dealing with the special case illustrated in Figure 2. In this case cell C is a coarse cell and either or both volume fluxes from C to its finer neighbours are positive. The fluxes then need to be computed independently for the top-and bottom-halves of the coarse cell (C t and C b in Figure 2). The equation of the interface in both halves can be determined easily using relations similar to ( 11), ( 12) and ( 13). Once this is done, calculation of the fluxes in each half-cell is identical to the standard case.

∆ t C a C u u i+1/2,j i+1/2,j
As has been noted before [START_REF] Rider | Reconstructing volume tracking[END_REF], this vof advection scheme is not strictly conservative because small over-or under-shoots of the volume fraction may occur for each direction-sweep.

Since interface reconstruction and geometrical flux calculation depend on consistent (i.e. bounded between zero and one) volume fractions, it is necessary to reset any volume fraction violating this constraint before proceeding to the next direction sweep, which will lead to loss of exact mass conservation. In practice we have found that this scheme is simple to implement on quad/octrees with mass conservation usually not being an issue, however it would be interesting to investigate the extension to quad/octrees of flux-redistribution schemes [START_REF] Rider | Reconstructing volume tracking[END_REF][START_REF] Lörstad | High-order surface tension VOF-model for 3D bubble flows with high density ratio[END_REF] or exactly mass-preserving geometrical schemes [START_REF] Scardovelli | Interface reconstruction with least-square fit and split eulerian-lagrangian advection[END_REF][START_REF] Aulisa | A geometrical area-preserving volume-of-fluid method[END_REF].

Conservation of interface spatial resolution

Depending on the criterion used for spatial adaptivity there may be cases when the interface spatial resolution could change during interface advection. Given two neighbouring cells C i and C i+1 this would happen whenever the following conditions are met:

1. C i+1 is coarser than C i , 2. u i+1/2 > 0, 3. C i+1 does not contain the interface (i.e. c i+1 = 0 or c i+1 = 1), 4. φ i+1/2 c i+1 .
In this case the interface fragment contained in cell C i will be advected to the coarser cell C i+1 resulting in a loss of spatial resolution of the interface.

In order to avoid this problem all cells verifying the conditions above are refined prior to interface advection for each of the direction sweeps. Initialising the volume fractions of the new children cells of C i+1 is trivial since by definition c i+1 = 0 or c i+1 = 1 (condition 3), however initialising the face-centred velocities u f of the new cells requires some thought. The face-centred velocities need to be discretely divergence-free (in order to ensure mass conservation during interface advection). They are usually obtained through the projection step (8) of section 2.2. In the present case however it would not be practical to reproject the whole face-centred velocity field just to get the new face-centred velocities in the few cells requiring refinement. A simple approach is to use a local projection to build the facecentred velocity field locally (Figure 3). In a first step the new face-centred velocities on the boundary of C i+1 (dashed or solid arrows) are initialised by copying the existing facecentred velocities of C i+1 and its neighbours (solid arrows). The second step initialises the remaining interior face-centred velocities u 0 , u 1 , v 0 and v 1 using an auxiliary potential p such that

u 0 ≡ p 1 -p 0 , u 1 ≡ p 3 -p 2 , v 0 ≡ p 0 -p 2 , v 1 ≡ p 1 -p 3 .
By construction the face-centred velocities in each fine cell will be conservative provided

u 0 -v 0 = d 0 , -u 0 -v 1 = d 1 , v 0 + u 1 = d 2 , v 1 -u 1 = d 3 ,
with d 0 , d 1 , d 2 and d 3 the sums of the incoming face-centred velocities on the boundaries of the corresponding fine cells (dashed or solid arrows). These conditions lead to the following discrete Poisson equation for field p

    -2 1 1 0 1 -2 0 1 1 0 -2 1 0 1 1 -2     p = d,
which has the solutions

p =     0 0 0 0 0 3/4 1/4 1/2 0 1/4 3/4 1/2 0 1/2 1/2 1     d + arbitrary constant,
provided the solvability condition d i = 0 is verified (i.e. the face-centred velocity field on the boundary of C i+1 must be divergence-free). This approach is easily generalised to three dimensions.

Balanced-force surface tension calculation

The accurate estimation of the surface tension term (σ κ δ s n) n+ 1 2 in the discretised momentum equation ( 3) has proven one of the most difficult aspect of the application of vof methods to surface-tension-driven flows.

In other methods such as front-tracking (e.g. the "pressure-gradient correction method" of [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF]) or levelset coupled with the "ghost fluid method" (gfm) [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF], the continuous description of the interface location can be used to derive accurate finite-volume or finite-difference estimates of the surface-tension and pressure-gradient terms taking into account the discontinuities at the interface. Unfortunately front-tracking methods cannot deal simply with topology changes and levelset methods have trouble conserving mass. Coupling vof, levelset and gfm can be used to overcome the individual drawbacks of each method [START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF][START_REF] Sussman | A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles[END_REF][START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF], however this introduces two different representations of the same interface with the associated complexity, efficiency and consistency issues.

In the context of vof methods, the original Continuum-Surface-Force (csf) approach of [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] proposes the following approximations

σ κ δ s n ≈ σ κ ∇c ˜, κ ≈ ∇ • n ˜with n ˜≡ ∇c |∇c ˜|,
where c ˜is a spatially-filtered volume fraction field. This approach is known to suffer from problematic parasitic currents when applied to the case of a stationary droplet in theoretical equilibrium [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF][START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF][START_REF] Harvie | An analysis of parasitic current generation in volume of fluid simulations[END_REF]. More recently Renardy & Renardy [START_REF] Renardy | PROST -A parabolic reconstruction of surface tension for the volume-offluid method[END_REF] and Francois et al. [START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF] noted that in this case, since the discretised momentum equation reduces to

-∇p n+ 1 2 + σ κ (δ s n) n+ 1 2 = 0,
Balanced-force surface tension calculation or equivalently using the csf approximation

-∇p n+ 1 2 + σ κ (∇c) n+ 1 2 = 0, (16)
it is possible to recover exact discrete equilibrium between surface tension and pressure gradient provided:

1. the discrete approximations of both gradient operators in equation ( 16) are compatible, 2. the estimated curvature κ is constant, then

p n+1/2 = σ κ c n+1/2 + arbitrary constant,
is the exact discrete equilibrium solution. An earlier implementation of a scheme with compatible pressure and interface gradient operators was also described in [START_REF] Sussman | An adaptive level set approach for incompressible two-phase flows[END_REF].

Condition 1 can easily be violated in naive implementations of the csf scheme [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF][START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF]. For example it may seem simpler to compute the surface tension force at the discrete locations where the volume fraction is defined. When using a classical staggered discretisation the surface-tension forces then need to be estimated on the faces of the control volume which is usually done by averaging the cell-centred estimates. The resulting discrete gradient operator applied to c in equation ( 16) is thus a spatially-averaged variant of the discrete gradient operator applied to p and condition 1 is violated.

In the present collocated scheme, the cell-centred pressure gradient is computed by averaging the face-centred pressure gradients as indicated in equation ( 9). The corresponding surface-tension steps verifying condition 1 are: a) apply the surface-tension force to the auxiliary face-centred velocity field u ⋆ f as

u ⋆ f u ⋆ f + ∆t σ κ n+ 1 2 f ρ(c n+ 1 2 f ) ∇ f c n+ 1 2 , (17)
b) apply the corresponding cell-centred surface-tension force to u ⋆ c as

u ⋆ c u ⋆ c + ∆t σ κ n+ 1 2 f ρ n+ 1 2 f ∇ f c n+ 1 2 c
.

These steps are applied immediately before the projection steps. In practice the implementation of steps a) and b) is identical to the implementation of the pressure corrections ( 8) and ( 9) respectively. Parameters are used to apply the routines either to c n+1/2 or p n+1/2 with weights σ κ n+1/2 or -1 respectively. This avoids duplicating similar code and guarantees that condition 1 is verified exactly.

Generalised height-function curvature calculation

Condition 2 depends on the accuracy of the curvature estimation. Estimating curvature has traditionally been the Achilles' heel of Volume-Of-Fluid schemes and has prompted many to recommend levelset [START_REF] Sussman | A level set approach for computing solutions to incompressible twophase flow[END_REF], coupled vof/levelset [START_REF] Sussman | A coupled level set and volume of fluid method for computing 3D and axisymmetric incompressible two-phase flows[END_REF][START_REF] Yang | An adaptive coupled level-set/volume-offluid interface capturing method for unstructured triangular grids[END_REF] or front-tracking schemes [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF][START_REF] Shin | Accurate representation of surface tension using the level contour reconstruction method[END_REF] as alternatives.

Besides the original curvature estimation of Brackbill, several methods have been proposed for estimating curvature from volume fraction fields. The height-function (hf) curvature calculation initially proposed by Torrey et al [START_REF] Torrey | NASA-VOF2D: a computer program for incompressible flows with free surfaces[END_REF] has recently been shown to give practical estimates of the curvature which are comparable in accuracy to estimates obtained using the differentiation of exact levelset functions [START_REF] Cummins | Estimating curvature from volume fractions[END_REF]. The hf method is also substantially simpler to implement than other high-order methods such as parabolic fitting (prost [START_REF] Renardy | PROST -A parabolic reconstruction of surface tension for the volume-offluid method[END_REF]) or spline interpolation [START_REF] Ginzburg | Two-phase flows on interface refined grids modeled with VOF, staggered finite volumes, and spline interpolants[END_REF]. An important shortcoming of the standard hf method which is rarely discussed, however, is that it becomes inconsistent when the radius of curvature of the interface becomes comparable to the mesh size. This case of an under-resolved interface often occurs in practice, particularly when dealing with topology changes.

In the following, I will present a height-function calculation generalised to quad/octree discretisations and which uses a hierarchy of consistent approximations to deal with the case of under-resolved interfaces.

Optimised height-function formulation on quad/octree discretisations

The standard height-function curvature calculation on 2d Cartesian grids [START_REF] Cummins | Estimating curvature from volume fractions[END_REF] can be summarised as:

1. Consider a 3 × 7 (resp. 7 × 3) stencil centred on the cell where the curvature is to be evaluated; an estimate of the interface orientation is used to choose the stencil best aligned with the normal direction to the interface. 2. Build a discrete approximation of the interface height y = h(x) (resp. x = h(y)) by summing the volume fractions in each column (resp. line). 3. Use finite-difference approximations of the derivatives of the discretised height-function to compute the curvature

κ = h ′′ (1 + h ′ 2 ) 3/2 x=0
While this scheme could easily be extended to a quad/octree discretisation using the approach described in section 3.1 (Algorithm 1), it would imply systematically reconstructing a large stencil around each cell where the curvature needs to be evaluated. In many cases the stencil does not need to be seven cells high (Figure 4.a) while some cases may require stencils higher than seven cells in order to obtain a consistent estimate of the interface height ( (a) A 3 × 3 symmetric stencil is sufficient to obtain a consistent estimate of the interface height for weakly-curved interfaces aligned with the grid. (b) Symmetric stencils higher than seven cells may be required to obtain consistent interface heights (here the 9 × 3 symmetric stencil indicated by the dotted lines is necessary). More compact asymmetric stencils can be built independently for each column (solid lines). (c) In this case only two of the three interface heights are consistent for either the vertical or the horizontal stencils.

An alternative to using wide symmetric stencils is to build local asymmetric stencils adapted to the geometry of the interface (in a manner similar to that proposed in [START_REF] Sussman | Improvements for calculating two-phase bubble and drop motion using an adaptive sharp interface method[END_REF]). For example in Figure 4.b local stencils of varying heights can be constructed independently for each of the three columns where the interface height is to be evaluated: a 1 × 6 stencil for the left column, a 1 × 3 stencil for the central column and a 1 × 4 stencil for the right column. The elementary operation is thus the evaluation of the interface height in a given column, for which I propose the following algorithm: C is a cell belonging to the column and can be seen as an initial guess of the position of the interface. The top and bottom directions used in steps 4.a and 8.a respectively are defined relative to the interface orientation so that for a consistent interface the cells at the top of the column are empty (c = 0) while the cells at the bottom are full (c = 1). In a first stage (steps 2 to 5) volume fractions in the top direction are summed until the interface has been found (I = true) and a full or empty cell has been reached. If the cell which has been reached is full a consistent stencil cannot be formed and an error is returned (condition 5). In a second stage (steps 6 to 9) the same procedure is repeated in the bottom direction, looking this time for a full cell. Given a good enough initial guess C this algorithm will be equivalent to using optimal local asymmetric stencils. Also, no a priori assumption is made on the maximum height of the stencil required to obtain a consistent interface height.

Algorithm 4 Interface height (C, top) 1. set H c(C) 2. set N C, c t c(C) 3. if c t < 1 then set I true else set I
Algorithm 4 is easily generalised to a quad/octree discretisation by considering virtual rather than actual cells in a manner similar to Algorithm 1. The virtual cells are constructed on demand when evaluating c(C) and c(N ).

For a given cell C and direction top, the full height-function curvature estimation algorithm can then be summarised as:

Interface height (N, top) c) set a common origin h

i h i + h(N i ) -h(N 0 )
3. if all heights are consistent a) return the curvature estimated using finite-difference approximations of the derivatives of the discretised height function h i

else a) return all the interface positions deduced from the consistent heights

Step 2.c is necessary since the stencils for each column are formed independently and do not necessarily share a common origin (cf. arrows in Figure 4.b). The function h(C) returns the absolute height of the centre of cell C.

Generalised height-function formulation for under-resolved interfaces

An important drawback of the standard height-function method is that even moderately curved interfaces can lead to configurations where consistent interface heights cannot be formed. An example of such a case is given in Figure 4.c where both the standard horizontal and vertical hf curvature estimations fail. The radius of curvature of the interface is 4 ∆ which is not a very coarse approximation for many practical applications.

While neither the horizontal nor the vertical stencils on their own can be used to construct a twice-differentiable discrete approximation of the interface height it is clear that a combination of all the stencils can allow such an approximation. In the case of Figure 4.c for example, four average interface positions can be obtained from the consistent interface height estimates (the circles in Figure 4.c). Fitting a curve (e.g. a parabola in 2d, paraboloid in 3d) through these points and differentiating the resulting analytical function will give an estimate of the curvature.

More generally given a cell C and n estimated interface positions {x 1 , , x n }, the following algorithm can be used to estimate the curvature

Algorithm 6 Parabola-fitted curvature (C, {x 1 , , x n }) 1.
if the number of independent interface positions is smaller than three (resp. six in 3d) a) a meaningful least-squares fit cannot be achieved, return an error 2. retrieve normal m to the interface in cell C (pre-computed using Algorithm 3) 3. compute the coordinates o of the barycentre of the reconstructed interface fragment contained in C 4. define an orthonormal coordinate system R ≡ {o, (i ′ , m)} (resp. {o, (i ′ , j ′ , m)} in 3d) 5. compute the transformed coordinates {x 1 ′ , , x n ′ } of the interface positions in R 6. fit a parabola (resp. paraboloid in 3d) by minimising

F (a i ) ≡ 1 j n z j ′ -f (a i , x j ′ ) 2 , with f (a i , x) ≡ a 0 x 2 + a 1 x + a 2 in 2d a 0 x 2 + a 1 y 2 + a 2 x y + a 3 x + a 4 y + a 5 in 3d 7. return the mean curvature at the origin o of R κ ≡          2 a 0 (1 + a 1 2 ) 3/2 in 2d 2 a 0 (1 + a 4 2 ) + a 1 (1 + a 3 2 ) -a 2 a 3 a 4 (1 + a 3 2 + a 4 2 ) 3/2 in 3d
The condition in step 1 uses the number of independent positions rather than the total number of positions n. Two interface positions a and b are considered independent provided |a -b| ∆ . This condition ensures that the resulting minimisation problem is wellconditioned. The minimisation only requires the solution of a 3 × 3 (resp. 6 × 6 in 3d) linear system.

In practice this scheme is sufficient to compute the curvature in most cases where the standard hf technique fails, however the number of independent positions -given by consistent interface height estimates -may become too small for complicated interfaces at very coarse resolutions. In this case a new set of interface positions is constructed by computing the coordinates of the barycentres of the reconstructed interface fragments in each cell of a 3 × 3 stencil (resp. 3 × 3 × 3 in 3d). If this approach in turn fails to provide enough independent positions, the cell considered most probably contains an isolated/degenerate interface fragment and the curvature is simply set to zero.

The overall algorithm for computing the interface mean curvature in a cell C containing the interface can then be summarised as The algorithm is hierarchical. In the majority of cases it will return within the first iteration of loop 3 (i.e. the standard hf method). In some cases the direction "best-aligned" with normal m will lead to inconsistent interface heights while the next best direction will not and the algorithm will return after more than one iteration of loop 3. Most other cases will be dealt with directly by step 5, rarer cases will require step 4.a and still rarer cases step 4.b. While estimating curvature using parabola fitting is substantially more expensive than the standard hf method, in practice the overall cost is thus largely dominated by the cost of the standard hf method.

The balanced-force surface tension correction [START_REF] Dai | Adaptive tetrahedral meshing in free-surface flow[END_REF] requires face-centred interface curvature estimates. The face-centred curvatures are computed either by averaging the cell-centred curvatures of the neighbouring cells when they both contain the interface, or by taking the value of the cell-centred curvature in either cell containing the interface.

Results

Convergence of the generalised height-function curvature

The standard height-function method has been shown to provide asymptotically second-order accurate estimates of the curvature when given an exact volume fraction field [START_REF] Cummins | Estimating curvature from volume fractions[END_REF]. In order to assess the consistency of the generalisation proposed in the previous section, we repeat the test case of a circular interface but over a larger range of spatial resolutions. The radius of the circle is varied from less than 2 ∆ to more than 16 ∆ and the corresponding volume fractions are initialised exactly. The position of the centre of the circle is varied randomly in the interval ([0, ∆], [0, ∆]), so that a representative sample of the various interface configurations is considered for each radius. The following relative error norms are defined

L 2 (κ exact ) ≡ 1 κ exact r i (κ r,i -κ exact ) 2 r i , L ∞ (κ exact ) ≡ 1 κ exact max r,i (|κ r,i -κ exact |),
where r is the index of the random circle position and i the index of a cell containing an interface fragment (i.e. for which 0 < c < 1). In addition two independent L ∞ norms are computed: one for the cells where the curvature is estimated using the standard hf method, the other for the cells using the parabola-fitting of Algorithm 6 (note that steps 4.a,b of Algorithm 7 are never necessary for this test case). The results obtained for 100 random circle positions (for each radius) are summarised in Figure 5. Second order First order Parabola % Fig. 5. Convergence of the relative curvature errors as a function of radius for a circular interface. Two sets of L ∞ (max) norms are given depending on whether the curvature was calculated using the standard height-function method (HF) or the parabola-fitting technique (parabola).

The thick curve in the figure indicates the proportion of cases where parabola fitting was necessary (from ≈ 80% for a radius of 2 ∆ down to ≈ 2% for a radius of 5 ∆). Above a radius of approximately 5 ∆, the standard hf method is always consistent and use of the parabola-fitting technique is not necessary. In this regime an asymptotically second-order convergence in both the L 2 and L ∞ norms is observed as in previous studies.

At low resolutions the largest errors are obtained for the cases in which the parabola-fitting technique is used. This is not surprising since at low resolution the hf method is consistent only for nearly horizontal or vertical interfaces for which good curvature estimates can be found. The parabola-fitting errors show second-order convergence. At intermediate resolutions (radius ≈ 4 ∆) the maximum errors for both methods are comparable and correspond to cases where the interface is inclined at 45 • (Figure 6). Interestingly, in this regime the L ∞ norm of the standard hf method only shows first-order convergence.

Over the whole range of resolutions these results are consistent with second-order convergence of the generalised height-function technique in both the L 2 and L ∞ norms.

Circular droplet in equilibrium 6.2.1 Stationary droplet

A circular interface with surface tension should remain at rest with the pressure jump at the interface exactly balancing the surface tension force (Laplace's law). In practice, depending on the method used for discretising the surface tension force and the pressure gradient, an exact numerical balance is difficult to obtain and often leads to so-called "spurious" or "parasitic" velocities fed by this imbalance [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF][START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF][START_REF] Harvie | An analysis of parasitic current generation in volume of fluid simulations[END_REF]. This severely limits the range of physical parameters which can be accurately simulated using csf-based techniques. The balancedforce surface tension formulation presented above is an important step toward an exact numerical solution of Laplace's problem, however an interesting recent study by Francois et al. [START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF] concluded that better methods for curvature estimation are required before this can be achieved. In contrast, in the following I will show that the combination of height-function curvature estimation and balanced-force surface tension described previously is sufficient to recover exact (to within round-off errors) numerical balance for Laplace's problem provided the shape of the interface is given enough time to relax to its numerical equilibrium shape. I will consider a circular interface centred on the top-left corner of the unit square. Symmetry conditions are applied on the top and left boundaries so that only a quarter of the droplet is simulated on a 32 × 32 grid. The diameter of the droplet is D ≡ 0.8 and the corresponding volume fractions are initialised exactly. The divergence criterion γ p is set to 10 -6 . The only velocity scale for the inviscid problem is

U σ ≡ σ ρ D ,
with ρ the constant density. This can be interpreted as the scale of the velocities associated with a capillary wave of wavelength comparable to D. The corresponding timescale is

T σ ≡ ρ D 3 σ ,
which is proportional to the period of the capillary wave. For a viscous fluid another timescale can be defined as

T ν ≡ D 2 /ν ,
with ν the kinematic viscosity. The time it takes for momentum to diffuse across the droplet is proportional to T ν . Finally the ratio of these two timescales is

T ν T σ = σ D ρ ν 2 = La √ ,
with La the Laplace number.

For our time-explicit discretisation of the surface tension term, numerical stability requires that the timestep be smaller than the period of the shortest capillary wave i.e. Figure 7 illustrates the evolution of the root-mean-square (rms) velocity with time for the range of Laplace numbers indicated in the legend. Time and velocity are made nondimensional using T ν and U σ as reference scales respectively. The rms velocity shows a roughly exponential decrease with time at all Laplace numbers. At long enough times (t ≈ T ν ) exact balance (to within round-off errors) is recovered between surface tension and pressure and the velocity tends toward zero. This result holds independently of the spatial resolution. Figure 8 gives a different representation of the initial oscillations apparent in Figure 7, with time made non-dimensional using T σ as reference scale. The case of an inviscid fluid (La = ∞) has also been added. The results for an inviscid fluid are comparable to those presented in Francois et al [START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF], who then concluded that this non-vanishing velocity was a measure of the spurious currents induced by an insufficiently accurate estimate of the curvature. I believe this is an incorrect interpretation of these results. From Figure 8 and looking at the details of the velocity field evolution (Figure 9 above and Figure 17.b,d in [START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF]) it is clear that the oscillations of the velocity field correspond to the oscillations of capillary waves of period ≈ 0.4 T σ (independent of the Laplace number). An estimate of the corresponding wavelength λ can be obtained using the relation [START_REF] Fyfe | Surface tension and viscosity with Lagrangian hydrodynamics on a triangular mesh[END_REF] T λ ≡ ρ λ 3 π σ ≈ 0.4 T σ , which gives

∆t ρ ∆ 3 π σ . (18)
λ ≈ 0.8 D ≈ π D/4,
which is consistent with the fourfold symmetry observed in the velocity field and in the initial curvature errors (case radius = 13 ∆ in Figure 6). These capillary waves exist because of the difference between the exact equilibrium interface shape (imposed as initial condition) and the numerical equilibrium shape (eventually reached at t ≈ T ν ). I suggest the correct interpretation is in fact that the velocity fluctuations observed consistent numerical solution of the evolution of a perturbed initial interface shape. The resulting capillary waves are then exponentially damped by viscosity -again in a physically consistent manner -on timescales of order T ν (with a damping coefficient proportional to ν as illustrated by the roughly constant slopes in Figure 7). This is in contrast with what is usually understood by spurious or parasitic currents which are continuously fed energy by an unphysical unbalance, do not vanish with time and cannot generally be interpreted as consistent physical entities (such as capillary waves).

When the numerical equilibrium interface shape is reached, we expect the estimated curvature to be exactly constant so that condition 2 of section 4 is verified. This is indeed the case as illustrated in Figure 10 where the saturation of the standard deviation at O (10 -11 ) is due to round-off errors. Finally the accuracy of the numerical equilibrium solution can be assessed by comparing the final shape of the interface with the exact equilibrium shape (the initial circular shape). Two measures are used: the relative curvature error |κ equilibriumκ exact |/κ exact and the shape error defined as

L 2 (shape) ≡ i (c i -c i exact ) 2 i , L ∞ (shape) ≡ max i (|c i
c exact is the exact volume fraction field. Figure 11 summarises the convergence of both measures when the spatial resolution is varied. Close to second-order convergence is obtained for both measures and the absolute values of the errors are small even at coarse resolutions. To conclude, the results presented in this section demonstrate that the combination of a balanced-force csf implementation and a height-function curvature estimation is sufficient to obtain an exact equilibrium solution (for velocity) in the case of a stationary droplet. The numerical equilibrium shape obtained is very close to the theoretical equilibrium shape and shows second-order convergence with spatial refinement. The timescale necessary to reach the numerical equilibrium solution is comparable to the viscous dissipation timescale T ν as expected from physical considerations. In practice this means that test cases designed to evaluate the accuracy of a given surface tension implementation (for the stationary droplet problem) need to make sure that simulations are run for timescales comparable to T ν (this was not the case in the study of Francois et al. [START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF]). For shorter timescales the results obtained will only reflect the accuracy of the initial curvature estimation. Estimating this initial error is better done using a test case similar to that of Section 6.1.

Translating droplet

While the previous test case demonstrates convergence toward the exact velocity solution in the case of a stationary droplet with surface tension, it is only a weak test of the combined accuracy of the interface advection scheme and the surface tension discretisation. In most practical applications droplets do not stand still but are advected around by the mean flow.

A simple but more representative variant can be built from the previous case by adding a constant, uniform background velocity field. The exact solution in the moving frame of reference is unchanged but the numerical solution requires advection of the interface across the fixed computational grid. Errors in interface advection estimates and hence -through surface tension -in the velocity field. Numerous studies have considered the case of a stationary droplet with surface tension [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF][START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF][START_REF] Renardy | PROST -A parabolic reconstruction of surface tension for the volume-offluid method[END_REF][START_REF] Shin | Accurate representation of surface tension using the level contour reconstruction method[END_REF][START_REF] Harvie | An analysis of parasitic current generation in volume of fluid simulations[END_REF][START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF] and/or the advection of a circular interface without coupling to the momentum equation [START_REF] Li | Calcul d'interface affine par morceaux[END_REF][START_REF] Rider | Reconstructing volume tracking[END_REF][START_REF] Scardovelli | Interface reconstruction with least-square fit and split eulerian-lagrangian advection[END_REF][START_REF] Yang | An adaptive coupled level-set/volume-offluid interface capturing method for unstructured triangular grids[END_REF][START_REF] Min | A second order accurate level set method on non-graded adaptive Cartesian grids[END_REF][START_REF] Malik | Adaptive VOF with curvature-based refinement[END_REF] but I am not aware of any study using a combination of both, hence the present test case.

The following results are obtained when considering a circular interface of diameter D = 0.4 contained in the unit square with a mesh refinement equivalent to a 64 × 64 regular spatial resolution. Periodic boundary conditions are used in the horizontal direction, symmetry conditions on the top and bottom boundaries. Figure 12 illustrates the evolution with time of the non-dimensional root-mean-square velocity -in the translating frame of reference -for different Laplace numbers. The Weber number is 0.4. Time and velocity are made non-dimensional using T U and U as reference scales respectively. In contrast to the results of the previous section, the velocity does not converge toward the exact solution with time. The frequency associated with the oscillations apparent in the figure scales like U /∆. This is due to advection errors continually perturbing the interface shape and thus the estimated curvature. The amplitude of this "spurious" velocity field is indicative of the accuracy of the combination of: advection scheme, curvature estimation and surface tension scheme. This simulation was repeated for a range of spatial resolutions and the convergence results are summarised in Figure ∞ spatial-norms are computed by taking the maximum of either norms over time. The maximum errors in velocity are of the order of 5% of U . Less than first-order convergence is observed for the maximum error and close to first-order convergence for the rms error. The shape errors can be computed using the exact translating interface shape as reference.

The maximum over time of the curvature and shape error norms are summarised in Figure 14. The roughly first-order convergence rates are consistent with the convergence rate of the error in velocity. Finally the dependences of the relative error in velocity on the Weber and Laplace numbers are illustrated in Figure 15.a and 15.b. The relative error scales like We -1/2 or equivalently like U -1 (σ being kept constant) which means that the absolute error in velocity is essentially independent from U . The relative error in velocity is also seen to be weakly dependent on the Laplace number (Figure 15 The convergence toward the exact solution for a stationary droplet presented in the previous section has been one of the main goals of several recent sophisticated surface tension implementations [START_REF] Renardy | PROST -A parabolic reconstruction of surface tension for the volume-offluid method[END_REF][START_REF] Shin | Accurate representation of surface tension using the level contour reconstruction method[END_REF][START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF]; however the results obtained in the present section -for a test case more relevant to real applications -underline the need for yet further advances. An avenue of investigation revealed by the translating droplet test case regards the tighter integration (and coupled evaluation) of the advection and curvature/surface tension schemes.

Capillary wave

Capillary waves are a basic feature of surface-tension-driven flows and their adequate numerical resolution is a prerequisite to more complex applications. The small-amplitude damped oscillations of a capillary wave are a now classical test case of the accuracy of numerical schemes for the time-dependent dynamics of viscous, surface-tension-driven twophase flows [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF]. For example a recent study by [START_REF] Gerlach | Comparison of volume-of-fluid methods for surface tension-dominant two-phase flows[END_REF] used this test case to assess the relative accuracies of a range of vof and vof/levelset techniques.

A sinusoidal perturbation is applied to a flat interface between two fluids initially at rest. Under the influence of surface tension the interface oscillates around its equilibrium position.

The amplitude of the oscillation decays due to viscous dissipation. The solution to this initial value problem is significantly different from the normal mode analysis of Lamb [START_REF] Lamb | Hydrodynamics[END_REF] (which predicts exponentially-damped oscillations) because the finite time of diffusion of vorticity generated at the interface into each of the phases introduces non-trivial history terms in the evolution equations. Nonetheless, Prosperetti [START_REF] Prosperetti | Motion of two superposed viscous fluids[END_REF] found an exact analytical solution in the limit of vanishingly small amplitudes which can be used as reference.

Prosperetti's theory is valid for infinite domains and we found that in order to get good agreement at high resolutions it is necessary to move the top and bottom boundaries far enough away from the interface. Using a [λ/2, λ/2][ -3 λ/2, 3 λ/2] domain -with λ the wavelength of the perturbation -seems to be sufficient. As in Popinet and Zaleski [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF] and Gerlach et al. [START_REF] Gerlach | Comparison of volume-of-fluid methods for surface tension-dominant two-phase flows[END_REF] the initial perturbation amplitude is λ/100, the Laplace number 3000 and the densities and viscosities of both fluids are identical. The timestep is controlled by the explicit surface-tension stability criterion [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF]. The typical evolution of the relative maximum interface elevation over time is illustrated in Figure 16. The normal-mode oscillation frequency ω 0 is defined by the dispersion relation

ω 0 2 = σ k 3 2 ρ ,
with k = 2 π/λ the wavenumber. The error between the theoretical solution of Prosperetti and the numerical solution can be evaluated using the L 2 norm

L 2 ≡ 1 λ ω 0 25 t=0 25/ω0 (h -h exact ) 2 ,
where h (resp. h exact ) is the maximum interface height obtained numerically (resp. theoretically). Table 1 and Figure 17 summarise the results obtained for: this study, the prost [START_REF] Renardy | PROST -A parabolic reconstruction of surface tension for the volume-offluid method[END_REF] and clsvof [START_REF] Sussman | A coupled level set and volume of fluid method for computing 3D and axisymmetric incompressible two-phase flows[END_REF] implementations of Gerlach et al [START_REF] Gerlach | Comparison of volume-of-fluid methods for surface tension-dominant two-phase flows[END_REF], the csf implementation of Gueyffier et al [START_REF] Gueyffier | Volume of fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF], the balanced-forced level-set method of Herrmann [START_REF] Herrmann | A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids[END_REF] and the fronttracking method of Popinet and Zaleski [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF]. Close to second-order convergence with spatial resolution is obtained for the present study, prost and clsvof. The deviation from second-order convergence for the csf and fronttracking methods at high resolutions may be due to the influence of the top and bottom boundaries (square domains were used in both these studies). The results of the present method are substantially more accurate at low resolutions than the results of the other methods given in Figure 17 and Table 1. This is important since in practice capillary waves are often the smallest features of interest and are usually captured at the limit of spatial resolution.

Method
At high resolutions the present method is at least as accurate as prost. The prost method relies on systematic parabola-fitting through non-linear minimisation which has been demonstrated to be very accurate [START_REF] Renardy | PROST -A parabolic reconstruction of surface tension for the volume-offluid method[END_REF] but is also much more computationally expensive than the present method.

Finally the convergence results for an air/water capillary wave are given in Table 2. The density ratio is 850 and the dynamic viscosity ratio is 55. [START_REF] Rayleigh | On the instability of a cylinder of viscous liquid under capillary force[END_REF]. The Laplace number in the denser fluid is 3000 and the other parameters are identical to the constant density case above. The density and viscosity are defined using the filtered field c ~of equations ( 1) and [START_REF] Li | Calcul d'interface affine par morceaux[END_REF]. Although the errors are larger than for the constant density case they are still small and good convergence is obtained. These results can be compared with the results in the recent study by Herrmann [START_REF] Herrmann | A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids[END_REF] (

Two-dimensional inviscid rising bubble

An inaccurate representation of surface tension is particularly noticeable for the case of a light gas bubble rising in a dense, low-viscosity fluid [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF][START_REF] Lörstad | High-order surface tension VOF-model for 3D bubble flows with high density ratio[END_REF][START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF]. Figure 18 illustrates the results obtained for one of the cases considered by Francois et al [START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF]. As is clear from the symbols in Figure their sensitive method pressure: circle symbols for the continuous approach (csf) and triangular symbols for the ghost-fluid sharp interface approach (ssf). The method described in the present article was applied at different resolutions with and without adaptive refinement. The solid line corresponds to the results obtained on the finest mesh using a constant resolution of 256 × 512. The dotted line gives the results for a constant resolution of 128 × 256 and the dashed line, the results for an adaptive resolution corresponding to an equivalent resolution of 256 × 512. The adaptive results were obtained by controlling the mesh size ∆ so that the vorticity criterion

∆ |∇ × u| < 1/2,
was verified at each timestep, with a limit on resolution of eight levels of refinement. A typical mesh is illustrated on Figure 18.(d). The vorticity generated at the interface is advected tangentially and accumulates at the stagnation points corresponding to the rear corners of the bubble. This increase in vorticity along the interface is matched by a corresponding increase in spatial resolution along the interface. The number and positions of the transitions between levels of refinement along the interface vary with bubble motion and shape evolution. This provides a good coverage of the consistency of the surface tension and curvature estimation schemes on the variable resolution quadtree mesh.

The results at t = 0.2 are identical for all methods and resolutions considered. At t = 0.35 the results for the current scheme are identical. They closely match the results of Francois et al in the front part of the bubble and fall between the csf and ssf results for the rear part. At t = 0.5, the results for the current scheme also fall between the csf and ssf results of Francois et al, but the coarse mesh results for the current scheme (dotted line) differ from the fine mesh results (a resolution dependence also observed by Francois et al). The constant fine mesh results and adaptive results are almost identical. This confirms the overall consistency of the adaptive scheme when resolution varies along the interface, for a flow with strong surface tension, low viscosity and large density ratio. Table 3 gives a summary of the mesh sizes and cpu times (on a single cpu 2.5 GHz 32 bits Intel Pentium processor) for the different cases considered. The timestep is controlled so that the explicit surface tension stability constraint [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF] is verified. The tolerance on the divergence of the velocity field γ p is set to 10 -3 . Setting a lower tolerance does not affect the results and mass is conserved to within less than 10 -2 %. A maximum of eight multigrid iterations is sufficient to reach this tolerance in all cases, one to three iterations are sufficient for 90% of the timesteps. The computational speed on the adaptive mesh is smaller than the speed on the regular mesh due mainly to the extra cost incurred by the more complex gradient operator at fine-coarse cell boundaries.

Three-dimensional capillary

is fundamental mechanism two-phase flow phenomena such as liquid/gas atomisation. It is also well-suited to dynamic adaptive mesh refinement as a wide range of flow scales need to be resolved, particularly around the breakup point [START_REF] Shi | A cascade of structure in a drop falling from a faucet[END_REF][START_REF] Brenner | Breakdown of scaling in droplet fission at high Reynolds number[END_REF][START_REF] Eggers | Nonlinear dynamics and breakup of free-surface flows[END_REF]. Dai and Schmidt [17] used their 3d adaptive unstructured tetrahedral mesh solver to study the case of the instability of a liquid cylinder with a free surface (i.e. the fluid motion of the gas phase was not resolved). Figure 19 reproduces their study (corresponding to Figure 16 of [START_REF] Dai | Adaptive tetrahedral meshing in free-surface flow[END_REF]) but for the case of a liquid jet surrounded by a lighter fluid. The computational domain is a unit box centred on the origin. Symmetry boundary conditions are applied for x = ± 1/2, y = -1/2, z = -1/2 so that only a quarter of the jet is simulated. Simple "outflow" boundary conditions are applied for y = 1/2 and z = 1/2 as

∂ n u n = 0, p = 0,
with n the direction normal to the boundary. The initial shape of the jet is given by radius

r(x) = r 0 [1 + ǫ sin(k x)],
with r 0 = 0.2, k = π and ǫ = 0.02. In order to compare with the solutions obtained from linear stability analysis, the velocity field is initialised as

u = ∇φ with φ ≡ B J 0 (i k r) sin(k x),
where the amplitude B is given by

≡ - ǫ r 0 c i k J 1 (i k r 0 )
, with c inviscid growth rate obtained by Rayleigh [START_REF] Rayleigh | On the instability of a cylinder of viscous liquid under capillary force[END_REF] as

c 2 ≡ σ ρ r 0 3 I 1 (k r 0 ) I 0 (k r 0 ) k r 0 (1 k r 0 2
).

The Laplace number ≡ σ r 0 /(ρ ν 2 ) is 238.34 as in [START_REF] Dai | Adaptive tetrahedral meshing in free-surface flow[END_REF]. The physical parameters for the surrounding fluid correspond to air/water density and viscosity ratios i.e. ρ g /ρ = 1.2/998 and µ g /µ = 1.8 × 10 -5 /1.003 × 10 -3 . The density and viscosity are defined using the filtered initial shape refinement levels. adjusted so that for any cell containing the interface the condition ∆ κ max < 1/5 is verified. The maximum level of refinement is set arbitrarily to ten. The timestep is continuously adjusted so that the stability condition ( 18) is verified.

Figure 20 illustrates the clean and continuous increase of resolution when approaching the breakup point. One can also note that -despite being essentially discontinuous -the vofreconstructed planar fragments provide a high-quality description of the interface. Fig. 20. Details of the interface just before ((t 0t)/T σ = 1.5 × 10 -3 ) and after breakup ((t 0t)/T σ = -8.5 × 10 -3 ). The interface is represented using the vof-reconstructed planar fragments in each cell. k r 0 = 0.628, ǫ = 0.02, La = 238.34.

The results are comparable to those of Dai and Schmidt although the present simulation is not limited to a singly-connected topology. As with any vof method no special treatment is required when breakup occurs and the simulation deals transparently with multiply connected interfaces (Figure 21). The maximum variation the total mass is less than 0.02%. Weber extended the study of Rayleigh to the case of a viscous fluid [START_REF] Weber | Disintegration of liquid jets[END_REF] and obtained the following equation for the growth rate of the perturbation

c 2 + c 3 ν k 2 - σ 2 ρ r 0 3 (1 -k 2 r 0 2 ) k 2 r 0 2 = 0.
Figure 22 summarises the results obtained for the evolution of the initial perturbation. Two numerical estimates of the perturbation are considered: the relative deformation associated with the maximum radial extent of the interface (r m axr 0 )/(r 0 ǫ), and the relative deformation associated with the minimum radial extent (r 0r m in )/(r 0 ǫ). Obviously the minimum radius vanishes at the time of breakup and the associated relative deformation takes a maximum value of 1/ǫ ( ≈ 50). The time is made non-dimensional using T σ ≡ (ρ r 0 3 /σ) For the maximum-radius relative deformation excellent agreement is obtained with either the inviscid theory of Rayleigh (La = 2000) or the viscous solution of Weber. Remarkably the growth rate is close to constant right up to the point of breakup. The good agreement with the inviscid theory at high Laplace number is also indicative of the overall low numerical dissipation of the discretisation scheme.

The minimum-radius relative deformation behaves similarly at initial times but departs from a constant growth rate before reaching the breakup point. The fluid motion asymptotically close to breakup can be shown to be a self-similar singularity of the Navier-Stokes equations with characteristic length scale

l ν = ν 2 ρ σ .
In this regime (r m in comparable to l ν ) the minimum radius is proportional to t 0t with t 0 the time of breakup [START_REF] Eggers | Universal pinching of 3D axisymmetric free-surface flow[END_REF][START_REF] Eggers | Nonlinear dynamics and breakup of free-surface flows[END_REF]. This regime is easily observed experimentally and numerically provided the ratio r 0 /l ν (i.e. the Laplace number) is small enough. Far enough away from breakup (l ν ≪ r m in ≪ r 0 ) and provided the Laplace number is large enough the fluid motion can be considered inviscid and scaling arguments give r m in ∼ (t 0t) 2/3 [START_REF] Peregrine | The bifurcation of liquid bridges[END_REF][START_REF] Eggers | Nonlinear dynamics and breakup of free-surface flows[END_REF].

Figure 23 shows the evolution of the minimum radius close to breakup for La = 2000 and 23.834. As expected two main regimes are clearly distinguished depending on the Laplace number. At La = 2000 the evolution of the minimum radius is accurately described by the inviscid theory for more than a decade of non-dimensional time. For La = 23.834 a transition between the two regimes occurs for r m in /r 0 ≈ l ν /r 0 = 1/La, the viscous regime (r min /r 0 < 1/La) being accurately described by the self-similar solution (r m in ∼ t 0t) while the inertial regime (r m in /r 0 > 1/La) is in reasonable agreement with the inviscid theory (r min ∼ (t 0t) 2/3 ).

La=2000

La=23.834 t 2/3 t ∆/r 0 Fig. 23. Evolution of the non-dimensional minimal radius r min /r 0 as a function of the non-dimensional time to breakup (t 0t)/T σ . The dash-dotted line (∆/r 0 ) indicates the limit of spatial resolution of the simulations.

The self-similar solution close to breakup predicts the formation of a thin liquid thread ("microthread" [START_REF] Kowalewski | On the separation of droplets from a liquid jet[END_REF]) connecting the two regions about to separate. Although these microthreads occur independently of the Laplace number, their characteristic length scale is l ν and they will be more easily seen at low Laplace numbers. Figure 24 illustrates the microthread obtained numerically just before breakup for La = 23.834. It looks very similar to the microthreads observed experimentally [START_REF] Kowalewski | On the separation of droplets from a liquid jet[END_REF][START_REF] Shi | A cascade of structure in a drop falling from a faucet[END_REF][START_REF] Eggers | Nonlinear dynamics and breakup of free-surface flows[END_REF]. The complex structure of the corresponding velocity field is given in Figure 24.b. Interestingly the solution has two stagnation points on the axis of symmetry (aside from the stagnation points on the symmetry planes at x = ± 1/2): one within the microthread (which will eventually lead to breakup) and a second one within the "neck" region. This result is consistent with the similarity solution obtained by Eggers [START_REF] Eggers | Universal pinching of 3D axisymmetric free-surface flow[END_REF] (see also Figure 17.b of [START_REF] Eggers | Nonlinear dynamics and breakup of free-surface flows[END_REF]). A recirculation region within the gas phase is associated with these two stagnation points. Figure 25 represents the evolution with time of the total number of leaf cells used to discretise the problem. The step increments in the number of cells (at t ≈ 7, 11, 12, ) correspond to the addition of one additional level of refinement when the curvature of the interface locally reaches the threshold ∆ κ max < 1/5. As expected the maximum number of cells ( ≈ 9,000) is reached at the time of breakup. After breakup the number of cells stabilises to about 2,500. This is larger than the initial number of cells ( ≈ 1,500) because the newlycreated, smaller droplets require a higher resolution (rightmost Figure 21). These numbers can be compared to the number of cells which would be required with:

1. a 3d static regular Cartesian grid with an equivalent resolution (ten levels of refinement):

(2 10 ) 3 ≈ 10 9 . 2. a 2d-axisymmetric static regular Cartesian grid: (2 10 ) 2 ≈ 10 6 . 3. a dynamically-refined octree grid with all the cells containing the interface refined at the same resolution ∆(t) controlled by ∆(t) max (κ m ax )(t) < 1/5: dashed curve in Figure 25.

It is clear that orders-of-magnitude savings in computation size can be achieved when using dynamic adaptation along the interface for this particular problem. Option 3 above which has been used in a number of recent studies [START_REF] Jeong | Finite element analysis of transient fluid flow with free surface using VOF (volumeof-fluid) method and adaptive grid[END_REF][START_REF] Greaves | A quadtree adaptive method for simulating fluid flows with moving interfaces[END_REF][START_REF] Kohno | Numerical analysis of moving interfaces a level set method coupled with adaptive mesh refinement[END_REF][START_REF] Wang | Finite-volume-type VOF method on dynamically adaptive quadtree grids[END_REF][START_REF] Theodorakakos | Simulation of sharp gas-liquid interface using VOF method and adaptive grid local refinement around the interface[END_REF][START_REF] Zheng | Adaptive unstructured remeshing-II: Application to two-and three-dimensional level-set simulations of multiphase flow[END_REF][START_REF] Greaves | Simulation of viscous water column collapse using adapting hierarchical grids[END_REF][START_REF] Yang | An adaptive coupled level-set/volume-offluid interface capturing method for unstructured triangular grids[END_REF][START_REF] Min | A second order accurate level set method on non-graded adaptive Cartesian grids[END_REF][START_REF] Nikolopoulos | Three-dimensional numerical investigation of a droplet impinging normally onto a wall film[END_REF] is simple to implement within a vof framework but does not make full use of adaptivity. The total cpu time ( ≈ 6,000 timesteps) for the simulation of Finally the three-dimensional results are cross-validated against axisymmetric results both with and without variable resolution along the interface. The axisymmetric results are obtained using a cylindrical-coordinates version of the code with k r 0 = 0.628, ǫ = 0.02, La = 2000. Two types of refinement are considered: a) A dynamically refined quadtree grid with all the cells containing the interface refined at the same resolution ∆(t) controlled by ∆(t) max (κ m ax (s))(t) < 1/5, with s the surface coordinate on the interface (equivalent to option 3 above). b) A dynamically refined quadtree grid with variable resolution along interface ∆(s, t) controlled by ∆(s, t) (s, t) < 1/5 (corresponding to the 3D simulation). Figure 26 summarises the results for these three cases: before breakup (a), to breakup and (d), and some time after breakup (c). The three-dimensional profile is obtained as a cross-section on the axis of of the 3D vof-reconstructed surface. The 3D solution is almost axisymmetric and the choice cross-section plane does not affect the results. The general agreement between the three solutions is excellent. Before breakup the axisymmetric results are identical. The small discrepancy between the 3D and axisymmetric results be due the fact the second radius of (in the plane perpendicular to the plane of axisymmetry) is computed completely differently in the axisymmetric and 3D cases. Breakup occurs slightly later the axisymmetric ∆(s, t) case which the somewhat larger discrepancy the axisymmetric results at later times (Figure 26.c). Table 4 gives the corresponding mesh sizes, number of timesteps and (per grid point and timestep). The for the 3d ∆(t) case are estimated the computing speed as in the 3d ∆(s, t) case. The axisymmetric ∆(s, t) and 3d cases save factors of 5 and 49 respectively, in mesh sizes and cpu time compared to the corresponding ∆(t) cases.

Discretisation

# grid points timesteps Speed 3d ∆(t) (estimated)

Conclusion

The combination of a quad/octree discretisation, csf surface tension scheme and height-function curvature estimation implemented the Gerris solver has been demonstrated to accurate and efficient solutions for surface-tension-driven flows. The method, height-function and csf formulations have been a fully-adaptive quad/octree refinement along interface. For the case of breakup of three-dimensional jet, leads to a reduction by factor of fifty of mesh compared to using a constant resolution along the interface. As a the dynamics of the breakup can be accurately and efficiently modelled across spatial range spanning the transition from inertial to regimes.

The height-function curvature technique has generalised to be consistent even at low resolutions. Test cases demonstrate second-order curvature across the whole range of spatial resolutions. Clean and robust surface-tension-driven interface breakup is demonstrated using this generalised technique.

In contrast to the recent study of Francois et al [START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF] the long-standing problem of "parasitic currents" around a stationary droplet is shown to be solved by the combination of appropriate implementations of a balanced-force csf approach and height-function curvature estimation. Exact balance (to machine accuracy) between surface-tension force and pressure gradient is recovered irrespective of spatial resolution and viscosity provided the interface is allowed enough time to relax to its equilibrium position. As expected this relaxation time is shown to be comparable to the viscous dissipation timescale. Furthermore the equilibrium shape of the interface shows second-order convergence towards the exact circular equilibrium shape as a function of spatial resolution.

Several issues require further investigation. The new "translating droplet in equilibrium" test case shows that the elimination of parasitic currents in the stationary case does not guarantee that the numerical solutions of other problems will be free from surface-tensioninduced velocity errors. Also, while the balanced-force csf method is shown to perform well it has the important shortcoming of not conserving momentum (in contrast to the continuum surface stress (css) approach of [START_REF] Gueyffier | Volume of fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF]). Initial tests suggest that this can be a major issue for problems involving long time-integrations and high surface-tension such as the rise of millimetric air bubbles in water.

Finally, I hope that providing the full implementation of the method as part of the freelymodifiable open source software Gerris [START_REF] Popinet | The Gerris Flow Solver[END_REF] will encourage other people to tackle these issues.

Fig. 1 .

 1 Fig. 1. Geometrical flux estimation.

Fig. 2 .

 2 Fig. 2. Special case for geometrical flux computation on a quadtree discretisation.

Fig. 3 .

 3 Fig. 3. Local construction of discretely divergence-free face-centred velocities.

  Figure 4.b).

  Fig. 4. (a)A 3 × 3 symmetric stencil is sufficient to obtain a consistent estimate of the interface height for weakly-curved interfaces aligned with the grid. (b) Symmetric stencils higher than seven cells may be required to obtain consistent interface heights (here the 9 × 3 symmetric stencil indicated by the dotted lines is necessary). More compact asymmetric stencils can be built independently for each column (solid lines). (c) In this case only two of the three interface heights are consistent for either the vertical or the horizontal stencils.

false 4 .

 4 while I = false or N contains the interface a) replace N with its top neighbour, update c t c(N ) b) update H H + c t c) if N contains the interface then set I true 5. if c t 0 then return inconsistent height 6. set N C, c b c(C) 7. if c b > 0 then set I true else set I false 8. while I = false or N contains the interface a) replace N with its bottom neighbour, update c b c(N ) b) update H H + c b c) if N contains the interface then set I true 9. if c b 1 then return inconsistent height 10. return H and N

Algorithm 5

 5 Height-function curvature (C, top) 1. set h 0 and N 0 Interface height (C, top) 2. foreach of the 2 (8 in 3d) columns neighbouring C a) define N as the closest neighbouring (virtual) cell to C in the plane perpendicular to the top direction b) set h i and N i

Algorithm 7

 7 curvature (C) 1. retrieve normal m to the interface in cell C (pre-computed using Algorithm 3) 2. compute a set S of two (resp. three in 3d) spatial directions in decreasing order of alignment with normal m 3. foreach top direction in S a) compute κ = Height-function curvature (C, top) b) if κ is consistent return κ c) else add the interface positions deduced from the consistent heights to set I 4. if the number of independent positions in I is smaller than three (resp. six in 3d) a) replace I with the set of interface positions built from the barycentres of the reconstructed interface fragments in a 3 × 3 (resp. 3 × 3 × 3 in 3d) stencil centred on Cb) if the number of independent positions in I is still too small return zero 5. return Parabola-fitted curvature (C, I)

Fig. 6 .

 6 Fig. 6. Distribution of relative curvature errors along the circular interface at two resolutions.

Fig. 7 .

 7 Fig. 7. Evolution of the rms velocity around a circular droplet in theoretical equilibrium for different Laplace numbers. Time and velocity are made nondimensional using T ν and U σ as reference scales respectively.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Initial evolution of the rms velocity around a circular droplet in theoretical equilibrium. Time and velocity are made non-dimensional using T σ and U σ as reference scales respectively.

Fig. 10 .

 10 Fig.10. Evolution of the standard deviation along the interface of the numerical curvature as a function of time for the Laplace numbers indicated in the legend. Time is made non-dimensional using T ν as reference scale.

Fig. 11 .

 11 Fig. 11. Convergence with spatial resolution of two error measures of the equilibrium interface shape. (a) Relative error on the equilibrium curvature: |κ equilibriumκ exact |/κ exact . (b) L 2 and L ∞ norms of the shape error.

Fig. 12 .

 12 Fig. 12. Evolution with time of the norm of the velocity in the translating frame of reference for different Laplace numbers. Time and velocity are made non-dimensional using T U and U as reference scales respectively.

Fig. 13 .

 13 Fig. 13. Convergence with spatial resolution of the non-dimensional velocity norms in the translating frame of reference. La = 12, 000.

Fig. 14 .

 14 Fig. 14. Convergence with spatial resolution of two error measures of the translating interface shape. (a) Relative error on the equilibrium curvature: |κ equilibriumκ exact |/κ exact . (b) L 2 and L ∞ norms of the shape error.

Fig. 15 .

 15 Fig. 15. Dependence of the relative error in velocity on: (a) the Weber number, La = 12, 000, (b) the Laplace number, We = 0.4 . Equivalent resolution 64 × 64.

Fig. 16 .

 16 Fig. 16. Time evolution of the relative maximum interface height of a capillary wave. La = 3000. Time is made non-dimensional using 1/ω 0 as reference timescale. Equivalent resolution 64 × 192.

Fig. 18 .

 18 Fig. 18. Shape of the bubble at times (a) 0.2, (b) 0.35 and (c) 0.5 for different methods. and triangles are the results obtained by Francois et al[START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF] (Figure21) using their csf and ghost-fluid (ssf) formulations respectively. The curve is the result of the current method for a resolution of 256 × 512, the dotted curve for constant resolution of 128 × 256 and the dashed curve for an adaptive resolution with eight levels of refinement. (d) Adaptive mesh and interface at t = 0.5.

Fig. 19 .

 19 Fig. 19. Initial stages of a capillary jet instability: k r 0 = 0.628, ǫ = 0.02, La = 238.34. Physical parameters correspond to a jet of water in air. The nondimensional times are from left to right 0, 11.32 and 12.1 (breakup time). The interface is represented using the vof-reconstructed planar fragments in each cell.

Fig. 21 .

 21 Fig. 21. Late stages of a capillary jet instability: k r 0 = 0.628, ǫ = 0.02, La = 238.34. The non-dimensional times are from left to right 12.31, 13.03 and 15.36.

Fig. 22 .

 22 Fig. 22. Evolution of the relative deformations of the interface. (a) (r maxr 0 )/(r 0 ǫ). (b) (r 0r min )/(r 0 ǫ).

Fig. 24 .

 24 Fig. 24. Detail of the interface close to breakup ((t 0t)/T σ = 1.5 × 10 -2 ) for a very viscous fluid: k r 0 = 0.628, ǫ = 0.02, La = 23.834. (a) vof-reconstructed interface. (b) Interface and streamlines in a cross-section on the axis of symmetry.

Fig. 25 .

 25 Fig. 25. Evolution of the total number of leaf cells required to discretise the jet while verifying the criterion ∆ κ max < 1/5. Maximum of ten levels of refinement. k r 0 = 0.628, ǫ = 0.02, La = 238.34.

Fig. 26 .

 26 Fig. 26. Comparison of axisymmetric and three-dimensional interface profiles. k r 0 = 0.628, ǫ = 0.02, La = 2000. t/T σ = 11.51. t/T σ = 11.61. (c) t/T σ = 11.79. Detail corresponding to the area in (b).

Table 1 .

 1 Convergence of the L 2 error on the maximum interface height of viscous capillary wave oscillations for different methods. The resolution is given in number of grid points per wavelength. ρ 1 /ρ 2 = 1, µ 1 /µ 2 = 1, La = 3000.

	/Resolution	8	16	32	64	128
	Present study	0.1568 0.0279 0.00838 0.0018 0.000545
	prost [66]	0.2960 0.0818 0.0069 0.0018	-
	clsvof [66]	0.3169 0.0991 0.0131 0.0033	-
	csf [19]	-	-	0.1168 0.0132 0.007
	rlsg [40]	-	0.1116 0.0295 0.0114 0.0067
	Front-tracking [15] 0.3018 0.0778 0.0131 0.0082 0.00645

Fig. 17. Convergence of the L 2 error on the maximum interface height of viscous capillary wave oscillations for different methods.

Table 15 )

 15 which show significantly larger errors and a saturation of convergence for resolutions above 64.

	Method/Resolution	8	16	32	64	128
	Present study	0.1971 0.0754 0.0159 0.00576 0.00313

Table 2 .

 2 Convergence of the L 2 error on the maximum interface height of viscous capillary wave oscillations for parameters corresponding to an air/water interface. ρ 1 /ρ 2 = 850, µ 1 /µ 2 = 55.72, La 1 = 3000.

Table 3 .

 3 Total number of grid points advanced, number of timesteps, cpu time (seconds) and computing speed (grid point × timestep / second) for the rising bubble problem on different discretisations.

  1/2 as reference timescale. Simulations were performed at three different Laplace numbers: 2000, 238.34 and 23.834.

	Relative deformation (maximum radius)	1 10 100	0	2	4	6	8	10 12 14 16 18 La=2000 La=238.34 La=23.834 Rayleigh Weber Breakup	Relative deformation (minimum radius)	1 5 10 20 50	0	2	4	6	8	10 12 14 16 18 La=2000 La=238.34 La=23.834 Rayleigh Weber Breakup
						Non-dimensional time						Non-dimensional time

Table 4 .

 4 Total number of grid points advanced, number of timesteps, cpu time (seconds) and computing speed (grid point × timestep / second) for the capillary breakup problem on different discretisations. k r 0 = 0.628, ǫ = 0.02, La = 2000.

Section 1

Section 4

Generalised height-function curvature calculation

 Section 5

 Section 6

 Section 6

Section 6

Acknowledgements

This work was funded by the New Zealand Foundation for Research, Science and Technology under contract no. CO1X0204.