Anton Smolianski
email: antsmol@amath.unizh.ch

Finite-Element/Level-Set/Operator-Splitting (FELSOS) Approach for Computing Two-Fluid Unsteady Flows with Free Moving Interfaces

Keywords: two-fluid interfacial flow, surface tension, Navier-Stokes equations, free moving boundary, finite element method, operator-splitting, level-set approach

Finite-element/level-set/operator-splitting (FELSOS) approach for computing two-fluid unsteady flows with free moving interfaces

Anton Smolianski

Introduction

Fluid flows with free moving surfaces or interfaces can be roughly divided into four general classes: bubbles/drops, jets, waves and films. Each class encompasses a large number of real-life physical phenomena having a great importance in diverse industrial applications. For example, bubble dynamics is of particular interest for chemical engineering, as bubbly flows are the core of bubble column chemical reactors; propulsion of liquid-metal jets constitutes the main part of metal forming processes; ocean waves are under thorough investigation in marine and coastal engineering, and liquid film flows are frequently encountered in coating and drying processes during paper or polymer production. It is worth noting that all above mentioned classes of fluid flows are, in essence, two-fluid flows, since even in the case when the second fluid is a gas (e.g., air) its dynamics cannot be neglected, with only few exceptions. Thus, in general, we have to deal with the flows of two immiscible fluids separated by their natural interface rather than with one-liquid free-surface flows.

There is a vast amount of literature devoted to numerical methods for free surface/interface fluid flows. As the comprehensive overviews containing a large number of references we would mention the papers [START_REF] Floryan | Numerical methods for viscous flows with moving boundaries[END_REF], [START_REF] Hou | Numerical solutions to free boundary problems[END_REF], [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF], and the book [START_REF] Shyy | Computational Fluid Dynamics with Moving Boundaries[END_REF].

Without loss of generality, the most popular way of classifying the numerical algorithms for fluid flows is to divide them into Eulerian, Lagrangian and mixed Eulerian-Lagrangian. Eulerian methods are characterized by a coordinate system that is stationary in the laboratory frame of reference. The fluid travels between different computational cells, in contrast to the Lagrangian methods, where each computational cell always contains the same fluid elements. Thus, Lagrangian methods are characterized by a coordinate system that moves with the fluid. The mixed Eulerian-Lagrangian methods rely on both Lagrangian and Eulerian concepts. This classification is very reasonable to describe the way of modeling the fluid flow, but does not contain any information on the approaches to modeling the interface motion. In this respect, there exists another commonly-used classification which treats all methods as either interface-tracking or interface-capturing. In the interface-tracking method the interface (free surface) is explicitly tracked along the trajectories of fluid particles in purely Lagrangian manner, which gives rise to the frequent use of interface-tracking in combination with Lagrangian or with mixed Eulerian-Lagrangian methods. The interface-capturing method is characterized by a reconstruction of the interface from the properties of appropriate field variables, e.g., fluid fraction or density. The latter classification clarifies the geometrical part of interfacial-flow modeling, i.e. the issues related to the interface motion. Finally, the way of coupling of the flow and the free moving interface can be "segregated" or "integrated": in segregated approach the flow is first computed with the "frozen" interface and, then, a new position of the interface is found using the last computed flow variables; in integrated approach the flow variables and new interface position are sought simultaneously.

To summarize, we may say that any computational method for free-surface flows consists of the following main ingredients:

(1) flow modeling (Eulerian, Lagrangian, mixed Eulerian-Lagrangian)

(2) interface modeling (tracking, capturing)

(3) flow-interface coupling (integrated, segregated)

The Lagrangian methods are naturally combined with the interface-tracking and the segregated flowinterface treatment. The strictly Lagrangian algorithm was used in [START_REF] Hirt | A Lagrangian method for calculating the dynamics of an incompressible fluid with free surface[END_REF] together with the finite-volume and in [START_REF] Okamoto | Two-dimensional sloshing analysis by Lagrangian finite element method[END_REF] with the finite-element method. In [START_REF] Fritts | The Lagrangian solution of the transient problems in hydrodynamics using a triangular mesh[END_REF] and [START_REF] Fyfe | Surface tension and viscosity with Lagrangian hydrodynamics on a triangular mesh[END_REF] the free Lagrangian method with remeshing/rezoning was employed, the meshless particulate Lagrangian approach was advocated in [START_REF] Monaghan | Simulating free surface flows with SPH[END_REF], [START_REF] Morris | Simulating surface tension with smoothed particle hydrodynamics[END_REF] (the Smoothed Particle Hydrodynamics (SPH) method) and the Boltzmann lattice-gas algorithm was addressed in [START_REF] Rothman | Lattice-Gas Cellular Automata[END_REF].

One of the most cited early papers within the Eulerian-Lagrangian framework is [START_REF] Hirt | An arbitrary Lagrangian-Eulerian computing method for all speeds[END_REF], where the algorithm called ALE (arbitrary-Lagrangian-Eulerian) was proposed. In this algorithm, the mesh is rearranged in the vicinity of the interface by letting the mesh nodes move with respect to the fluid in a prescribed manner; away from the interface the grid is kept fixed. The advective velocity is correspondingly corrected to take into account the grid motion, and the interface is tracked by following the Lagrangian motion of vertices aligned initially with the interface. The flexibility in dealing with the motion of mesh vertices makes the ALE-type methods very attractive for free-surface flow simulations, and the algorithms of this type were successfully used in [START_REF] Hughes | Lagrangian-Eulerian finite element formulation for incompressible viscous flow[END_REF], [START_REF] Belytschko | Finite element methods with user-controlled meshes for fluid-structure interaction[END_REF], [START_REF] Donea | An arbitrary Lagrangian-Eulerian finite element method for transient dynamic fluid-structure interaction[END_REF], [START_REF] Ramaswamy | Arbitrary Lagrangian-Eulerian finite element method for the analysis of free surface fluid flows[END_REF], [START_REF] Maury | Characteristics ALE method for unsteady free surface flows with surface tension[END_REF], [START_REF] Bänsch | Numerical Methods for the Instationary Navier-Stokes Equations with a Free Capillary Surface[END_REF]. All these works relied on the finite-element method and on the segregated treatment of flow-interface coupling. The ALE methodology was exploited also in [START_REF] Hansbo | The characteristic streamline diffusion method for the time-dependent incompressible Navier-Stokes equations[END_REF] and in [START_REF] Tezduyar | A new strategy for finite element computations involving moving boundaries and interfaces -The deforming-spatial-domain/space-time procedure[END_REF], where the space-time finite-element method was combined with least-squares type stabilization.

There is a special group of methods based on Lagrangian-Eulerian conception of mesh movement and on the fully coupled ("integrated") treatment for the system "flow variables -interface". In these methods the system of flow equations with free-surface boundary conditions is discretized as a whole with respect to the flow variables and to some functional representation (parametrization) of the interface. The resulting system of the nonlinear algebraic equations is then solved using a Newton or quasi-Newton iterative procedure. Such methods were proposed in [START_REF] Ruschak | A method for incorporating free boundaries with surface tension in finite element fluid-flow simulators[END_REF] and [START_REF] Saito | Study of coating flow by the finite element method[END_REF] for steady free-surface flows, and then extended to unsteady flows with free moving boundaries in [START_REF] Engelman | Finite element simulation of incompressible fluid flows with a free/moving surface[END_REF], [START_REF] Cuvelier | A time dependent free boundary governed by the Navier-Stokes equations[END_REF], [START_REF] Christodoulou | Discretization of free surface flows and other moving boundary problems[END_REF].

Eulerian methods are used in combination with either interface-tracking or interface-capturing approach. The former approach can be further decomposed into surface-tracking and volume-tracking. Surface tracking methods represent an interface as a series of interpolated curves through a discrete set of points on the interface. At each time step, the information about the location of the points and sequence in which they are connected is saved. The points are then moved according to an interface evolution equation. For the overview of early works on surface-tracking methods, the paper [START_REF] Hyman | Numerical methods for tracking interfaces[END_REF] may be consulted. The later works using surface-tracking approach are due to Glimm's group (see [START_REF] Glimm | Front tracking applied to Rayleigh-Taylor instability[END_REF], [START_REF] Glimm | Robust computational algorithms for dynamic interface tracking in three dimensions[END_REF]) and Tryggvason's group [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF], [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF]). Among the recent works on surface-tracking we could mention [START_REF] Shyy | Computational Fluid Dynamics with Moving Boundaries[END_REF], [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF] and [START_REF] Tornberg | Interface Tracking Methods with Application to Multiphase Flows[END_REF]. The modified surface-tracking methods not using a connectivity of the set of interface points have been proposed in [START_REF] Torres | The point-set method: front-tracking without connectivity[END_REF] and in [START_REF] Shin | Modeling three-dimensional multiphase flow using a level contour reconstruction method for front tracking without connectivity[END_REF]. The comparison of front-tracking with lattice-Boltzmann method has been recently carried out in [START_REF] Sankaranarayanan | A comparative study of lattice Boltzmann and front-tracking finite-difference methods for bubble simulations[END_REF]. In all these works, tracking of the interface was computationally segregated from the calculation of flow variables.

Volume-tracking Eulerian methods do not store a representation of the interface but reconstruct it whenever necessary. The reconstruction is done cell by cell and is based on the presence of marker quantity within the cell. The marker particles are only used to show which cells contain fluid and are moved with a fluid velocity in a purely Lagrangian manner. The first Eulerian volume-tracking algorithm for free-surface flows seems to be the marker-and-cell (MAC) method of [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF] that used fixed uniform mesh and the finite-difference approximation. The approach was extended and strenghtened by many researchers; we would mention here the papers [START_REF] Hirt | Calculating three-dimensional flows around structures and over rough terrain[END_REF] and [START_REF] Ramshaw | A numerical technique for low-speed homogeneous two-phase flow with sharp interfaces[END_REF]. For more recent algorithms using marker particles idea the works Glowinski et al. (1992) and [START_REF] Nakayama | An Eulerian finite element method for time-dependent free surface problems in hydrodynamics[END_REF] may be consulted, in which the finite element method was employed for the spatial discretization. It is worth noting that, like in other methods using purely Eulerian way of flow modeling, the treatment of flow-interface coupling is performed in a segregated manner within the volume-tracking framework.

In the interface-capturing Eulerian methods the interface is reconstructed from the properties of suitable field variables, such as fluid fractions. Namely, the interface is represented as either a discontinuity line of some characteristic function ("discontinuous interface-capturing") or a zero-level set of some implicit function ("continuous interface-capturing"). That function obeys pure transport equation, which states that the interface is a material line propagating with the fluid.

The first algorithm of the discontinuous interface-capturing type was suggested in [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] and is called volume-of-fluid (VOF) method. The VOF-type algorithms typically employ a segregated treatment for the system "flow variables -interface" and finite difference or finite volume approximations on fixed grids. For the review on state-of-the-art VOF-like methods one may be pointed to the papers [START_REF] Rudman | Volume-tracking methods for interfacial flow calculations[END_REF] and [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF] (see also [START_REF] Scardovelli | Interface reconstruction with least-square fit and split Eulerian-Lagrangian advection[END_REF]). Some further references include the work [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] remarkable by its "continuum surface force" (CSF) approach proposed to include the surface tension into the right-hand side of momentum equation (see also [START_REF] Liovic | Numerical modelling of free surface flows in metallurgical vessels[END_REF]) and the paper [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF] where the alternative "continuum surface stress" method was proposed within the finite volume framework (see also [START_REF] Wu | Simulation of two-fluid flows by the least-squares finite element method using a continuum surface tension model[END_REF] for the implementation of the latter method in the combination with finite elements). The VOF interface-capturing method has been exploited in the absence of surface tension effect also in [START_REF] Kawarada | Numerical method for a free surface flow on the basis of the fictitious domain method[END_REF], [START_REF] Kelecy | The development of a free surface capturing approach for multidimensional free surface flows in closed containers[END_REF], [START_REF] Puckett | A high-order projection method for tracking fluid interfaces in variable density incompressible flows[END_REF], [START_REF] Aliabadi | Stabilized-finite-element/interface-capturing technique for parallel computation of unsteady flows with interfaces[END_REF].

In contrast to the representation of the interface as a discontinuity line within discontinuous interfacecapturing framework, in the continuous approach the interface is defined as a zero level set of some continuous function. This gives rise to the notion of "level-set approach" that started from the work [START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF] and has been further developed in [START_REF] Sussman | A level set approach for computing solutions to incompressible twophase flow[END_REF], [START_REF] Chang | A level set formulation of Eulerian interface capturing methods for incompressible fluid flows[END_REF], [START_REF] Tornberg | Interface Tracking Methods with Application to Multiphase Flows[END_REF], [START_REF] Quecedo | Application of the level set method to the finite element solution of two-phase flows[END_REF], [START_REF] Codina | A numerical model to track two-fluid interfaces based on a stabilized finite element method and the level set technique[END_REF], [START_REF] Iafrati | A domain decomposition approach to compute wave breaking (wave-breaking flows)[END_REF], see also the books [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF] and [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF]. A very similar "pseudo-concentration" technique was developed in a parallel manner, see [START_REF] Thompson | Use of pseudo-concentrations to follow creeping viscous flows during transient analysis[END_REF], [START_REF] Codina | Mould filling simulation using finite elements[END_REF], [START_REF] Lewis | Efficient mould filling simulation in castings by an explicit finite element method[END_REF], [START_REF] Lock | Local mesh adaptation technique for front tracking problems[END_REF].

The detailed comparison of diverse numerical methods for interfacial flows led us to some specific choice of basic components of the numerical modeling strategy. First, we choose purely Eulerian approach, since it enables us to use a fixed structured grid on a fixed computational domain. Second, we rely on the interfacecapturing in order to be able to deal with complex interfacial motions including interface merger, folding and break-up. In particular, the level-set approach is taken in the present work. Next, we employ the operatorsplitting approach that immediately yields a segregated treatment of not only the flow-interface coupling but also of different parts of the problem, which correspond to different physical processes. This is computationally very advantageous, as the specialized numerical scheme can be used for each part of the problem's operator, and, instead of one very large problem, we have to resolve a sequence of smaller subproblems.

The rest of the paper is organized as follows. In Section 2 the mathematical model is set out, Section 3 is devoted to the description of the computational method and the discussion of corresponding stability issues, Section 4 presents diverse numerical experiments, and, finally, the conclusions are drawn in Section 5.

Setting of the problem

We start with the specification of the physical assumptions used to derive a suitable mathematical model for our problem. We consider an unsteady laminar flow of two immiscible fluids. Both fluids are assumed to be viscous and Newtonian. Moreover, we suppose that the flow is isothermal, thus neglecting the viscosity and density variations due to changes of a temperature field. We assume also that the fluids are incompressible. The validity of this assumption is affected by several factors, the most important of which is the condition for Mach number to be smaller than, approximately, 1/3 (see Batchelor (1967), Ch. 3.6, for a thorough discussion on the incompressibility assumption). That condition is satisfied in our case, since we deal with essentially subsonic flows. Presuming, in addition, the fluids to be homogeneous, we may infer that the densities and viscosities are constant within each fluid. We utilize the sharp-interface (zero interfacial thickness) approach; the density and viscosity have, therefore, a jump discontinuity at the interface (see, e.g., [START_REF] Batchelor | An Introduction to Fluid Dynamics[END_REF]). We assume that the interface has a surface tension. We also suppose that there is no mass transfer through the interface (i.e. the interface is impermeable), and there are no surfactants present in the fluids (hence, there is no species transport along the interface). Under such conditions we do not have to consider the variations of surface tension coefficient in tangential to the interface direction, i.e. the solutocapillary Marangoni effect (the thermocapillary Marangoni effect has been excluded by the assumption on isothermal character of the flow). Therefore, the surface tension coefficient may be assumed constant.

Suppose that the motion of two viscous immiscible fluids under our investigation is confined to some box (a parallelepiped in 3D, a rectangle in 2D). For the sake of simplicity we consider only 2D case, but the proposed approach is easily extendible to 3D situation. The boundary of the box can be physical (e.g., the walls of a container), artificial (if we consider a flow in unbounded domain) or partly artificial. In fact, we can always restrict ourselves to some bounded region of interest and consider, then, the flow in that region only. On the other hand, this enables us to avoid dealing with asymptotics at infinity and to make the problem more tractable from computational viewpoint. We denote the boundary of the box by Σ, the domains occupied with the fluids by Ω 1 and Ω 2 and the interface between the fluids by Γ (Γ = ∂Ω 1 ∩ ∂Ω 2 , where ∂Ω i is the boundary of Ω i , i = 1, 2), see Figure 1. Let also Ω be the entire region occupied with the fluids, i.e. the interior of the box (Ω = Ω 1 ∪ Ω 2 ∪ Γ). The domains Ω 1 and Ω 2 may be multiply connected, and the interface Γ may intersect the box boundary Σ. Taking into account the physical assumptions considered above, we may assert that the flow of both fluids is governed by the incompressible Navier-Stokes equations

ρ(x) ∂v ∂t + v • ∇v -∇ • (2µ(x)S) + ∇p = ρ(x)g , (1)
∇ • v = 0 in Ω , t > 0 (2)
supplemented with the interfacial conditions

[v] Γ = 0 , [-pI + 2µS] Γ • n = κσn . (3)
Here the density ρ(x) = ρ 1 in Ω 1 and ρ 2 in Ω 2 , the viscosity µ(x) = µ 1 in Ω 1 and µ 2 in Ω 2 , g is the acceleration of gravitational field, S = 1 2 (∇v + (∇v) T ) is the deformation rate tensor, I is the identity tensor, σ is the coefficient of surface tension, κ is twice the mean curvature of the interface, n is the unit normal to the interface, and [...] Γ denotes a jump across the interface Γ.

The equations ( 1)-( 2) and interfacial conditions (3) should be complemented with some boundary condition on Σ for velocity, for example, v = 0 on Σ (4) and with the initial conditions

Γ t=0 = Γ (0) (5) v t=0 = v (0) in Ω (0) , (6)
where Γ (0) is the initial position of the interface determining initial shapes of the domains Ω (0)

1 and Ω (0)

2 . Obviously, we have three unknowns here: the velocity of fluid v(x, t), the pressure p(x, t) and the interface Γ(x, t). The position of the interface at any moment of time can be determined using the fact that, in the absence of mass transfer through the interface (see the first of the interfacial conditions (3)), the interface is convected by the fluid, i.e. the interface normal velocity is equal to the normal component of fluid velocity (the latter is continuous across the interface, see (3)).

The problem at hand contains three key ingredients: (i) flow equations (i.e. the Navier-Stokes equations with discontinuous coefficients and singular capillary force), (ii) moving interface and (iii) coupling between velocity-pressure fields and the interface (through the coefficients, capillary force and interfacial advective velocity). To attack the problem numerically we advocate the operator-splitting approach; namely, at each time step, we first resolve the Navier-Stokes system with fixed known interface, then, using computed velocity field, we find the new approximation of the interface. Having found the new interface position we can calculate its normal and curvature, and, thus, evaluate the surface tension force and the density/viscosity coefficients to be used on the next time step in the Navier-Stokes equations.

3 Numerical method

Discretization of the Navier-Stokes equations

Using the Marchuk-Yanenko fractional-step scheme (see, e.g., [START_REF] Marchuk | Methods of Numerical Mathematics[END_REF]) we may separate the convective nonlinearity, viscous diffusion and incompressibility from one another and treat each of them with corresponding numerical technique. Thus, on each time interval [t n ; t n+1 ] the Navier-Stokes (NS) system (1)-( 3) is approximated by a sequence of three subproblems:

1. NS-convection step ∂v ∂t + v • ∇v = 0 in Ω × (t n ; t n+1 ) , (7)
v t=tn = v n , =⇒ v * . 2. Viscous step ρ(x) ∂v ∂t -∇ • (2µ(x)S) = ρ(x)g in Ω × (t n ; t n+1 ) , (8)
[v] Γ = 0 , [2µS] Γ • n = κσn on Γ , v t=tn = v * , =⇒ v * * . 3. Projection step ∂v ∂t + 1 ρ(x) ∇p = 0 , (9)
∇ • v = 0 in Ω × (t n ; t n+1 ) , [v] Γ • n = 0 , [p] Γ = 0 on Γ , v t=tn = v * * , =⇒ v n+1 , p n+1 .
Here v n (x) and p n (x) are the approximations to v(x, t n ) and to p(x, t n ) , n = 0, 1, . . ., respectively, and v 0 ≡ v (0) . Each of three subproblems must be complemented by suitable boundary conditions which will be discussed below. The general strategy consists in resolving the steps successively; moreover, each step can be treated with its own time discretization scheme employing a specific (variable) time-step size dictated by stability and/or accuracy reasons. Thus, the given global time step ∆t = t n+1 -t n may be subdivided into smaller time steps within each of the subproblems.

The fractional-step scheme is formally 1st-order accurate but known to be very robust (see, e.g., [START_REF] Glowinski | Finite element methods for Navier-Stokes equations[END_REF]). A rather similar operator-splitting approach was proposed in Zienkiewicz & Codina (1995), Zienkiewicz et al. (1995) (Characteristic-Based Split (CBS) method) and proved to be very efficient for simulation of transient fluid flows.

The operator-splitting implies a decomposition of the interfacial conditions. The splitting of the interfacial stress jump condition indicates that the surface tension balances the jump of viscous stress only (subproblem 2), while the pressure is simply continuous across the interface (subproblem 3). Such approximation of the interfacial condition seems to be in contradiction with the fact that the pressure does have a jump discontinuity at the interface, which can be clearly seen from the famous Laplace-Young equation, if the viscosity is very small. However, we will demonstrate in numerical tests that chosen scheme is capable of capturing the pressure interfacial discontinuity with a good accuracy; thus, there is no need to enforce the discontinuity explicitly. It is worth noting that such kind of splitting of the interfacial stress condition is common for operator-splitting (projection) methods applied to the Navier-Stokes equations with interfaces (see, e.g., [START_REF] Sussman | A level set approach for computing solutions to incompressible twophase flow[END_REF] for the finite-difference or [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF] for the finite-volume implementation).

All the subproblems ( 7)-( 9) are discretized in space by the finite element method using a single uniform triangulation for both velocity and pressure. It is worth to emphasize that the grid is kept fixed throughout the whole process of computation, independently of the interface location. If we denote the triangulation by T h , where h is the mesh size (i.e. the maximal diameter of the triangles), the approximation space for the pressure can be defined as

Q h = {p h ∈ C 0 (Ω) p h T ∈ 1 ∀T ∈ T h } , ( 10
)
where 1 is the space of polynomials in two variables of degree less than or equal to 1. The discrete velocity belongs to the following space:

V h = {v h ∈ C 0 (Ω) v h T ∈ * × * ∀T ∈ T h } , (11)
where * is defined as (a

) * ≡ 1 or (b) * ≡ b 1 . Here b 1 = 1 span{λ 1 λ 2 λ 3 } (λ i , i = 1, 2, 3
, are the area (barycentric) coordinates on the triangle T ). In the case (a) we deal with the 1 / 1 finite element for the velocity and pressure, and in the case (b) with the so-called "mini" finite element (see, e.g., [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF]). It is worth noting that the "mini" element satisfies the LBB (inf-sup) condition, while the equal-order ( 1 / 1 ) approximation does not. However, in combination with the operator-splitting ( 8)-( 9) (which is, in fact, a continuous version of the well-known Chorin projection scheme), the 1 / 1 approximation is stable, if the time step ∆t ≥ Ch 2 (C is a constant), see, e.g., [START_REF] Guermond | On stability and convergence of projection methods based on pressure Poisson equation[END_REF]. This fact enables us to use the most economical but sufficiently accurate 1 / 1 approximation for the velocity and pressure. Now we address briefly the solution of the subproblems ( 7)-( 9). The Navier-Stokes convection step is the first-order nonlinear hyperbolic problem, and, thus, the boundary conditions consist of prescribing the fluid velocity only on the inflow part of the boundary Σ (i.e. where the fluid velocity is directed inward the domain Ω). Then the problem is solved with the explicit Taylor-Galerkin scheme TTG-4A of [START_REF] Quartapelle | High-order Taylor-Galerkin methods for nonlinear multidimensional problems[END_REF].

The viscous step requires the boundary conditions to be as for the complete problem, i.e. ( 4). The use of the finite element method for spatial discretization of the viscous step is very essential, since in the weak (variational) formulation, on which the method relies, the interfacial stress jump condition becomes a natural condition. Thus, after multiplying (8) with weighting function w and integrating in Ω by parts, the surface tension force becomes automatically incorporated into the variational formulation:

Ω ρ(x) ∂v ∂t • w dx + Ω 2µ(x)S • •∇w T dx = Ω ρ(x)g • w dx + Γ κσn • w dΓ . ( 12
)
The method does not require to approximate a delta-function as it is usually done for accounting of singular capillary force (see, e.g., [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF], [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], [START_REF] Sussman | A level set approach for computing solutions to incompressible twophase flow[END_REF]) and, hence, precludes the interface smearing. As we are about to see, the interface Γ can be easily localized owing to the finite element discretization; thus, a sharp interface can be maintained, and the capillary force can be computed by the direct integration over Γ. Another advantage of the considered approach consists in alleviating the singularity connected with the differentiation of discontinuous viscosity coefficient; in fact, this problem is completely eliminated with the formulation (12). After an implicit temporal discretization of ( 12) with the backward Euler or the Crank-Nicolson scheme (the term "implicit" is used for the velocity field only, while the interface Γ and, consequently, the density and viscosity coefficients are treated explicitly), we obtain a linear algebraic system having symmetric and positivedefinite matrix; the system could be successfully solved, for instance, by incomplete Cholesky preconditioned conjugate gradient method.

The projection step defines an inviscid flow problem, thus, the boundary conditions consist in prescribing the boundary values only for the normal component of fluid velocity on Σ. After implicit time-discretization with the time step ∆t the problem can be recast in the form of a Poisson-type equation for the pressure endowed with the homogeneous Neumann boundary condition

-∇ • ( 1 ρ(x) ∇p n+1 ) = - 1 ∆t ∇ • v * * in Ω , (13)
n • ∇p n+1 Σ = 0 , (14)
where v * * is the intermediate velocity obtained on the viscous step. The final divergence-free velocity v n+1 is then derived as

v n+1 = v * * -∆t 1 ρ(x) ∇p n+1 . ( 15
)

Approximation of the interface

We use the level-set approach (see [START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF], [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF]) for capturing the interface Γ. The approach consists in specifying a continuous "level-set function" Φ such that Φ is the signed distance to the interface Γ, Φ > 0 in Ω 2 and Φ < 0 in Ω 1 . Obviously, we have Γ = {x| Φ(x, t) = 0 ∀t ≥ 0} and |∇Φ| = 1 in Ω.

Since the interface is convected with the fluid velocity, the following evolution equation for Φ(x, t) can be easily derived (see, e.g., [START_REF] Chang | A level set formulation of Eulerian interface capturing methods for incompressible fluid flows[END_REF]):

∂Φ ∂t + v • ∇Φ = 0 in Ω , t > 0 . ( 16
)
We solve this pure advection equation with the fluid velocity v found as the solution of the Navier-Stokes system and with the same explicit Taylor-Galerkin scheme TTG-4A of [START_REF] Quartapelle | High-order Taylor-Galerkin methods for nonlinear multidimensional problems[END_REF] as on the Navier-Stokes convection step. An explicit scheme with the CFL stability condition Co ≤ 1 (Co = |v|∆t/h is the Courant number), particularly the scheme TTG-4A, is well suited for solving the level-set convection problem, as this condition means that during one time step the interface does not move farther than to the distance h from its previous position. Thus, the interface does not skip any domain element, which is proved to be a reliable guide in selecting the time step for interface propagation (see, e.g., [START_REF] Lock | Local mesh adaptation technique for front tracking problems[END_REF]).

The continuous piecewise-linear approximation is utilized for Φ on the same uniform grid that is used for velocity-pressure approximation. Such approximation allows us to easily localize the interface at any moment of time: indeed, the intersection of the approximate level-set function with each triangle of the grid is just a line segment possibly degenerated to a point (a vertex of the triangle); moreover, to find this line segment one has to determine only the intersection points of the approximate, linear level-set function with triangle's edges. Thus, using piecewise linear continuous approximation for the level-set function on simplicial mesh, we obtain unique piecewise linear representation for the interface Γ.

Level-set reinitialization

At the beginning of computational process the level-set function Φ is initialized as the signed distance to the interface Γ (0) . But, as time goes, the level-set function may become very steep or flat in some regions, particularly in the vicinity of the interface, thus losing the nice property |∇Φ| = 1. This deterioration of the level-set function is a natural consequence of the convection process, but it makes difficult an accurate determination of the interface. In order to cure the situation, the level-set function must be reinitialized, i.e. made again the signed distance function. It is worthwhile to note that this "redistancing" should be done only in some vicinity of the interface, since the values of the level-set function far from Γ have no influence on the interface dynamics. Using the interface localization discussed above, the redistancing can be done by straightforward computation of the normal distance to the interface at each grid point located within the mh-band around the interface (integer m ranges typically between 5 and 10). The algorithm has a complexity O(N 2 ), where N is the number of grid nodes along one coordinate direction, i.e. depends linearly on the total number of grid points (see [START_REF] Smolianski | Numerical Modeling of Two-Fluid Interfacial Flows[END_REF]).

Level-set correction

The reinitialization of the level-set function is very important for accurate determination of the interface, but it cannot guarantee the mass conservation. In fact, due to the incompressibility assumption the area (volume in 3D) of the region occupied by each of the two fluids must be conserved during the whole computational process. However, this is not the case because of different numerical errors, and after many time-steps or any change of the interface topology the relative decrease/increase of the area of one of the fluids may be about several percents (see, e.g., [START_REF] Sussman | A level set approach for computing solutions to incompressible twophase flow[END_REF]). We propose a very simple and efficient approach that uses an additional, level-set correction step to explicitly enforce the mass conservation. Suppose that we have performed the reinitialization of the level-set function Φ, thus, in some neighbourhood of the interface, Φ is the signed distance function and its level sets (isolines) are equidistant. Then we may easily correct the area of the domain Ω 2 = {x ∈ Ω : Φ > 0}, still retaining the shape of the interface, by considering as a new zero-level set some isoline of Φ that lies in the vicinity of the current interface curve (see Figure 2). This operation can be accomplished simply by moving the level-set function upward or downward, i.e. by adding to Φ some signed constant C Φ , where |C Φ | is the distance between old and new zero-level sets (the fact, that Φ is a distance function, plays the important role here).

Φ = Φ + new Γ Γ c c new new Φ Φ Ω Ω 2 2
Figure 2: Level-set correction (solid lines are the old and the new zero-level sets).

To find the expression for C Φ , we denote by Φ new the new level-set function, Ω new 2 = {x ∈ Ω : Φ new > 0}, and utilize the well-known formula for the first variation of a volume integral (see, e.g., [START_REF] Delfour | Shape derivatives and differentiability of Min Max[END_REF] or [START_REF] Cuvelier | Some numerical methods for the computation of capillary free boundaries governed by the Navier-Stokes equations[END_REF]):

S exact -S(Ω 2 ) = Ω new 2 dx - Ω2 dx = Γ (C Φ n) • n dΓ + O(C 2 Φ ) = C Φ Γ dΓ + O(C 2 Φ ) , (17)
where S exact is the exact area of the region occupied with the second fluid (S exact is always known to us), S(Ω 2 ) is the area of Ω 2 . It immediately follows that

C Φ = S exact -S(Ω 2 ) L(Γ) (18)
where L(Γ) is the length of the interface Γ. The formula ( 18) is accurate up to O(C 2 Φ ) (we assume that 0 < l 1 ≤ L(Γ) ≤ l 2 < +∞, l 1 and l 2 are independent of any physical and numerical parameters of the problem).

First, it is noteworthy that if S(Ω 2 ) > S exact we obtain from (18) C Φ < 0, and the level-set function Φ is to be lowered, which automatically implies a shrinkage of the domain Ω 2 ; if S(Ω 2 ) < S exact then C Φ > 0, and Ω 2 automatically expands. Next, we may note that the level-set correction procedure makes sense also in the case of multiply-connected domain Ω 2 and in the case when Γ is only a part of the boundary of Ω 2 .

The justification of the level-set correction is essentially based on two facts: Φ being a distance function in some vicinity of the interface and smallness of |C Φ |. The former fact means that the level-set correction should follow the reinitialization on every time-step. The smallness of modulus of the correction constant is the crucial point in justifying the method; indeed, artificial moving of the interface is not allowed, except as the movement magnitude is not greater than the interface interpolation error. If we assume that Γ ∈ C 2 , the piecewise linear interfacial interpolation gives us the error O(h 2 ) in L ∞ -norm. Hence, |C Φ | should not be greater than O(h 2 ); below we will show that it can be expected. It is worth noting that the effect of the level-set correction is usually almost negligible within one time-step, but the correction's main purpose is to prevent an accumulation of numerical errors in a long run.

The major sources of the error in the interface location are the discretization error of the convection scheme and the discretization error of the advective velocity. The Taylor-Galerkin scheme TTG-4A employed for the level-set convection is 4th-order accurate in time and has the 4th-order spatial accuracy at grid nodes (see Thomee & Wendroff (1974) and [START_REF] Donea | An introduction to finite element methods for transient advection problems[END_REF]). Since we use 1 -elements for the level-set function Φ h , which immediately yields a piecewise linear representation for Γ h , the L ∞ -norm of the error in the interface convection should not be greater than O(h 2 ), if the exact interface Γ is smooth enough (the temporal error does not pollute the spatial accuracy, since ∆t = O(h) due to accuracy and stability reasons). The error coming from the discrete advective velocity is also not greater than O(h 2 ): in fact, the interface propagates with the fluid velocity, whose discretization error is not worse than O(h) in L ∞ -norm (as will be shown in the static bubble test, see also [START_REF] Smolianski | Numerical Modeling of Two-Fluid Interfacial Flows[END_REF]), thus, during one time-step of size ∆t = O(h) the error in the interface location cannot become greater than O(h) • O(h) = O(h 2 ). Hence, we can expect the magnitude of zero-level set movement during each mass-correction to be not greater than O(h 2 ).

Remark 3.1. Since the error of formula ( 18

) is O(C 2 Φ ) and |C Φ | is O(h 2
), the level-set correction method should imply the mass conservation up to O(h 4 ). However, in reality we deal with the function Φ h , which is not an exact distance function but only its piecewise linear approximation; thus, moving Φ h upward or downward in accordance with C Φ we obtain a new zero-level set that may differ by O(h 2 ) from the desired one. This implies the mass conservation error O(h 2 ) which is consistent with the spatial accuracy of the interface approximation.

Approximation of interface normal and curvature

Having computed the piecewise-linear approximation Φ h of the level-set function, we can find the interface normal by virtue of the formula n h = ∇Φ h /|∇Φ h |. This results in a piecewise-constant approximation of the normal; if we assume that the interface Γ (hence, the function Φ in some neighbourhood of the interface) belongs to C 3 , a standard gradient averaging procedure applied to n h will lead to the piecewise-linear continuous approximation ñh having O(h 2 ) accuracy (see [START_REF] Křižek | Superconvergence phenomenon in the finite element method arising from averaging gradients[END_REF], [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF], [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique[END_REF]). It is important that, since Φ h is defined over the entire domain Ω, ñh (x) is also formally defined for almost every x ∈ Ω. Using this fact, we may compute the curvature by two methods: either directly as

κ h = -∇ • ñh (19)
or by using the variational formulation for the curvature equation

Ω κh q h dx = Ω ñh • ∇q h dx - Σ (ñ h • n Σ )q h dΣ , κh ∈ Q h , ∀q h ∈ Q h , (20)
where n Σ is the outward unit normal to the boundary Σ and Q h is the finite-dimensional space of continuous piecewise-linear functions defined on the considered uniform triangulation (see ( 10)). Solution of (20) simply amounts to resolving the system with the consistent mass matrix. It is demonstrated in [START_REF] Smolianski | Numerical Modeling of Two-Fluid Interfacial Flows[END_REF] that the piecewise-constant approximation κ h has O(h) accuracy, while the piecewise-linear κh attains the accuracy O(h 2 ). The latter fact is quite remarkable, since we use only linear interface approximation, and allows us to evaluate accurately the surface tension force.

Evaluating the interfacial force and density/viscosity fields

After the approximations for the interface normal and curvature have been constructed, the discrete interfacial force term

Γ h κh σ ñh • w h dΓ h (21)
can be easily computed (κ h stands here for either κ h or κh ). Since Γ h is the union of closed interfacial line segments {γ hi }, the integral ( 21) can be represented as a sum of the integrals over γ hi and, thus, computed elementwise. The integral over each interfacial line segment is evaluated with the help of the standard 2-point Gaussian quadrature that is exact for polynomials of the 3rd degree. This choice of the quadrature formula is justified if the weighting function w h for the viscous momentum equation is piecewise linear ( 1 -elements for velocity). Numerical experiments show that the same quadrature rule can be utilized when b 1 -elements are used for velocity, although the integration will not be exact.

The localization of the interface enables us to compute accurately the integrals containing discontinuous density or viscosity coefficients; this is especially important for preserving a discrete momentum balance in the vicinity of the interface. Define the characteristic function

I(x) = 1 for x ∈ Ω 1 0 for x ∈ Ω 2 ,
and the discontinuous density and viscosity fields

ρ(x) = ρ 2 + I(x)(ρ 1 -ρ 2 ) , µ(x) = µ 2 + I(x)(µ 1 -µ 2 ) .
It is evident that on the elements not crossed by the interface, the density and viscosity coefficients are simply constant, and evaluation of the corresponding integrals can be done in a standard way. Thus, the problem is to compute the integrals of the form

Ti ρ(x)A(x) dx (22)
over each triangle T i intersected by the interface Γ h (A is a continuous within T i function). First, we may note that taking ρ(x) constant in T i (equal, e.g., to the average (ρ 1 + ρ 2 )/2) implies the error O(h 2 ) of computing the integral over T i , and, since the number of such triangles crossed by Γ h is proportional to 1/h, the error in computing the global integral Ω ρ(x)A(x) dx is O(h) for a general function A. This is insufficient accuracy, and numerical experiments indicate that this method is approximately equivalent to the smoothing out the density/viscosity coefficients over a few grid cells in the finite difference framework (see the discussion in [START_REF] Tornberg | Interface Tracking Methods with Application to Multiphase Flows[END_REF]). On the other hand, a straightforward use of standard numerical quadratures for a discontinuous integrand is strictly prohibited, as the quadratures are not convergent in such a case.

To overcome the difficulty of integrating a discontinuous function, we advocate the so-called "discontinuous integration" approach that is proved to be useful in two-phase problems (see, e.g., [START_REF] Fachinotti | A fast convergent and accurate temperature model for phasechange heat conduction[END_REF] and references herein, [START_REF] Tornberg | Interface Tracking Methods with Application to Multiphase Flows[END_REF]). The approach essentially uses the fact that the discrete interface Γ h is just a line segment within each interfacial element T i . Indeed, the triangle T i becomes cut in two parts, a triangle T * i and a quadrilateral. If we denote by I * and I * * the values of the characteristic function I(x) in those parts of T i , the integral ( 22) can be calculated as follows:

Ti ρ(x)A(x) dx = (ρ 2 + I * (ρ 1 -ρ 2 )) T * i A(x) dx +(ρ 2 + I * * (ρ 1 -ρ 2 )) Ti\T * i A(x) dx = (I * -I * * )(ρ 1 -ρ 2 ) T * i A(x) dx + (ρ 2 + I * * (ρ 1 -ρ 2 )) Ti A(x) dx .
This expression can be readily computed, since I = 1 where the level-set function Φ < 0 and I = 0 where Φ > 0 (we set Φ positive in Ω 2 and negative in Ω 1 ), and both integrals are over triangles; no subdivision of the quadrilateral T i \T * i into triangles is required, which makes the algorithm simpler (this is especially true for its 3D counterpart). On the discrete level the integrals over T i and T * i may be computed exactly, assuring a good evaluation of the discrete balance equations near the interface.

The error of the discontinuous integration method comes from the interface approximation by a straight line, which yields O(h 3 ) error of the integration over a single triangle T i (with an error constant depending on A(x)). The number of the triangles intersected by the interface is proportional to 1/h, thus the error in computing the global integral Ω ρ(x)A(x) dx is O(h 2 ) for a general function A.

Remark 3.2. It is shown in [START_REF] Ferziger | Interfacial transfer in Tryggvason's method[END_REF] that numerical smoothing (smearing) of the physical parameters like density and viscosity across the interface may lead to significant errors in the balance equations. The presented combination of the finite element method and the level-set approach allows to avoid any numerical smearing of material coefficients or surface tension force at the interface.

Summary of the algorithm

Our computational approach for numerical modeling of interfacial flows can be summarized as follows:

Step 0. Initialization of the level-set function and velocity.

For each n-th time-step, n = 1, 2, . . .:

1. Computing of interface normal, curvature and density/viscosity fields. 2. Navier-Stokes convection step. 3. Viscous diffusion step. 4. Projection step. 5. Level-set convection step. 6. Reinitialization step. 7. Level-set correction step.

The steps 1.-7. are performed successively, and each of the steps 2.-5. may use its own local time-increment size. On each step the last computed velocity is exploited; the viscous diffusion and projection steps use the interface position found on the previous global time-step. It is also noteworthy that the steps 5.-7. of the n-th time-step can be computed in a fully parallel manner with the step 2. of the next, (n + 1)-st time-step. The whole algorithm is very flexible; it allows, for instance, to compute unsteady interfacial Stokes flow just by omiting the Navier-Stokes convection step.

Stability issues and time scales

Each interfacial flow problem is characterized by some set of non-dimensional parameters (criteria) as well as by some specific time scale. That time scale is extremely important, as it reflects the dynamics of the considered physical process; it also serves as a guideline for selecting a time-step size in numerical simulations. Although the physical time scale is always dependent on the problem at hand, there are some common scales connected with typical physical phenomena like gravity, surface tension, viscosity.

Let L be the length scale of the considered problem. The gravity phenomenon is well represented by the dimensionless Froude number F r = U (g) √ gL , where g is the magnitude of the gravitational acceleration and U (g) is the characteristic scale of velocity. If we denote by ∆t (g) phys the typical time scale of a gravity-driven physical process, then U (g) = L/∆t , where µ is the dynamic viscosity, σ is the coefficient of surface tension and U (ca) is the velocity scale. If ∆t (ca) phys is the typical time scale of a capillarity-driven physical process, then U (ca) = L/∆t Finally, the viscosity is known to give a "viscous time-scale" through the non-dimensional Reynolds number

Re = ρL 2 µ∆t (v) phys , if L/∆t (v)
phys is considered as the velocity scale. Requiring Re ≈ 1 we immediately obtain ∆t

(v) phys ≈ ρL 2 µ .
To capture the dynamics of the modelled physical process, the numerical method should use a time-step not exceeding the time-scales of all involved physical phenomena. However, besides physical restrictions on the time-step size, there are some constraints originated from the numerical approximation. In our computational method, an example of such stability-restriction is the CFL condition appearing on the Navier-Stokes convection and on the level-set convection steps due to explicit time-discretization. The projection step does not suffer from any upper bound for the time increment, having only the inverse stability-constraint (lower bound) for the time-step size if the equal-order interpolation is used for velocity and pressure (this constraint is excluded by using an LBB-stable approximation). The viscous diffusion step possesses two stability constraints due to explicit treatment of the gravity and surface tension forces. Those constraints can be easily derived by the following empirical analysis (see [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF]):

the stability restriction due to gravity can be represented in a form of CFL-type condition

vg ∆t (g) num h
≤ 1, where the corresponding gravity-induced velocity v g = g∆t (g) num ; then, we obtain the maximal allowed time-step size as ∆t (g) num = h g ;

the stability restriction due to capillarity can be represented in a form of CFL-type condition

vca∆t (ca) num h ≤ 1
, where the corresponding capillarity-induced velocity v ca can be calculated using the capillarity-induced acceleration δ h σκ h ρ (δ h is the discrete delta-function included since the capillary force appears only on the interface); it is easy to see that δ h = 1 h at the interface, and the maximal possible curvature resolved on the grid with size h is κ h = 1 h ; thus, we obtain

v ca = σ h 2 ρ ∆t (ca)
num , and finally

∆t (ca) num = ρ σ h 3/2 .
The latter condition was derived also in [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF].

If an implicit time-discretization for the viscosity term is used, there are no other stability restrictions except those indicated above. Thus, the final stability-constraint for the viscous diffusion step is

∆t v ≤ min{∆t (ca) phys ; ∆t (v)
phys ; ∆t (g) num ; ∆t (ca) num } , since ∆t

(g)
phys is always greater than ∆t (g) num .

Remark 3.3. It is interesting to note that ∆t (g) num and ∆t

(ca)
num observed in numerical experiments are almost 2 times larger than their theoretical predictions, which can be explained by the stabilizing influence of the level-set correction procedure (see Section 3.2).

Numerical results

The static bubble test

In order to verify our numerical approach, we start with the simple test considered by [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF] (see also [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF]): a circular bubble in static equilibrium. In this problem, the net surface force should be zero, since at each point on the bubble surface the tension force is counteracted by an equal and opposite force at a diametrically opposed point. The correct solution is a zero velocity field and a pressure field that rises from a constant value of p out outside the bubble to a value of p in = p out + σ R inside the bubble, according to the Laplace-Young law (R is the bubble radius).

Assume that the computational domain is the square (0, 1) × (0, 1); the bubble has the center (0.5; 0.5) and the radius R = 0.25. We set the viscosity µ and density ρ equal to 1 everywhere in the domain and impose the no-slip condition on the domain boundary; the pressure is fixed by setting it zero at one of the domain corners (hence, p out = 0 in our computations). The gravity is neglected; the surface tension coefficient σ will be varied, thus yielding different values of the Laplace number La = (2R)σρ/µ 2 . The computations are performed until the non-dimensional time t = t/∆t (ca) phys = 250, where the physical (capillary) time scale ∆t

(ca) phys = (2R)µ/σ.
As a result of numerical simulations one can observe spurious velocity currents of amplitude U (see Figure 3). In the work by [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with SURFER[END_REF] it was conjectured that the amplitude of the spurious currents must be proportional to σ/µ, which is equivalent to having an approximately constant value of the capillary number Ca = U µ/σ. The conjecture on constancy of the capillary number is clearly confirmed in our computations, see Table 1 illustrating the constant character of Ca over a broad range of Laplace numbers. Table 1: Constancy of the capillary number Ca with respect to the Laplace number La (mesh size h = 1/40, 1 elements for velocity).

Table 2 shows the convergence of the non-dimensional velocity in the l ∞ norm, i.e. the decrease rate for the non-dimensional amplitude U µ/σ of the spurious currents (the l ∞ norm is defined as the absolute maximum over the grid nodes). The experiments were performed on three different grids with h = 1/20, 1/40, 1/80; the Laplace number was fixed at 5.0 • 10 3 . Two types of elements were tested for the velocity approximation: 1 andb 1 .

1 b 1 h = 1/20 1.4 • 10 -2 1.5 • 10 -2 h = 1/40 9.1 • 10 -3 9.2 • 10 -3 h = 1/80 5.0 • 10 -3 4.2 • 10 -3 p ≈ 0.8 0.9
Table 2: The error and the order of convergence (p) in the l ∞ norm for the non-dimensional velocity v h µ/σ: comparing the two types of velocity elements.

It is apparent that both 1 and b 1 elements exhibit nearly the 1st order convergence. The velocity resulting h = 1/20 6.9 • 10 -4 8.6 • 10 -4 h = 1/40 1.8 • 10 -4 2.3 • 10 -4 h = 1/80 4.7 • 10 -5 5.9 • 10 -5 p ≈ 1.9 1.9

Table 3: The error and the order of convergence (p) in the discrete L 1 norm for the non-dimensional velocity v h µ/σ: comparing the two types of velocity elements.

In Table 3 the convergence of the non-dimensional velocity in the discrete L 1 ("average") norm is shown (the discrete L 1 norm is defined as

||v h || 1 = 1 M M i=1 |v h (i)|,
where {v h (i)}, i = 1, M is the set of nodal values of the discrete velocity, and M = O(h -2 ) is the total number of velocity nodes). Since all the errors are concentrated within an interface neighbourhood covering O(1/h) nodes, the rate of convergence in the discrete L 1 norm should be one order higher than in the l ∞ norm. This is confirmed by the numerical results of Table 3. Figure 4 illustrates the pressure behavior in dependence on the mesh size. The 1 / 1 finite element pair was used for the velocity-pressure approximation; the b 1 / 1 elements gave qualitatively the same results. It is clear that, in spite of the continuous approximation for the pressure, the pressure discontinuity at the interface is maintained well. We have to note that the proposed operator-splitting algorithm is better suited for purely transient problems, while the steady-state results may suffer from the splitting error (which is of order of ∆t = O(h), see, e.g., [START_REF] Glowinski | Finite element methods for Navier-Stokes equations[END_REF]). The truly sharp discontinuity of the pressure will be observed below (see Figure 12) in the numerical experiments with a rising bubble. The qualitatively good agreement with the Laplace-Young law is supported by the quantitative results shown in Table 4. The absolute and relative errors in the Laplace-Young law indicate very good accuracy of the pressure inside of the bubble (p in ) when the 1 elements are employed for velocity approximation; the b 1 velocity elements give nearly the same pressure accuracy.

|p in -p out -σ R | |pin-pout-σ R | (σ/R) • 100% h = 1/20 2.0 • 10 -2 0.5% h = 1/40 2.0 • 10 -3 0.05% h = 1/80 6.0 • 10 -4 0.015%
Table 4: The absolute and relative errors in the Laplace-Young law.

To summarize the results of the presented test, we may say that the splitting of the interfacial stress jump condition implies, as it could be expected, some error in the velocity gradient, but the error is concentrated in a small neighbourhood of the interface. While the velocity approaches zero like O(h) in the l ∞ norm, the velocity gradient remains of order O(1) in this norm to balance the surface tension force at the interface. In the L 1 norm the convergence becomes one order higher, which clearly indicates the local character of the spatial distribution of the errors. The velocity interfacial perturbation resulting from the operator-splitting makes the pressure gradient at the interface negative in the radial direction (i.e. away from the bubble's center); this, in combination with the zero velocity (and, thus, pressure gradient) inside and outside of the bubble, gives rise to the correct pressure profile, see Figures 3 and4. It is worth noting that the absence of any oscillations near the pressure discontinuity line can be attributed to the stabilizing character of the projection-step operator: the pressure is sought as the solution of the Poisson equation that, in its turn, acts as a "natural" diffusion operator at the interface.

Rising bubble

The next test is the simulation of a single bubble, rising in an initially quiescent fluid due to the effects of buoyancy. The problem is characterized by four dimensionless parameters, which we choose to be the density ratio ρ 1 /ρ 2 , the viscosity ratio µ 1 /µ 2 , the Reynolds number Re = (2R) 3/2 √ gρ 1 /µ 1 and the Eötvös number Eo = 4ρ 1 gR 2 /σ. Here R is the initial radius of the bubble, index "1" corresponds to the fluid surrounding the bubble, index "2" to the fluid inside of the bubble. The Morton number M = Eo 3 /Re 4 is often used instead of the Reynolds number.

The computational domain is the rectangle (0; 1) × (0; 2); the bubble is initially circular, with the center (0.5; 0.5) and the radius R = 0.25. The free-slip boundary condition is imposed on the vertical walls, and the no-slip condition on the horizontal walls; the initial velocity is zero everywhere in the domain.

Figures 5 and6 show the evolution of the bubble shape and of the velocity field, corresponding to the same density and viscosity ratios (ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 ). The Reynolds number Re = 35 in both cases, but the Eötvös number equals 25 and 125 respectively. The computations were performed on the triangulation based on 40 × 80 rectangular grid, the time-step size ∆t = h/2. The Crank-Nicolson time discretization was employed for viscous diffusion term, and the 1 elements were used for the velocity.

It is clearly seen that when the surface tension coefficient is rather large (Eo = 25) the bubble has an ellipsoidal shape, while smaller surface tension (Eo = 125) results in a more significant deformation of the bubble. In the latter case the bubble has a shape of a "dimpled ellipsoidal-cap". Although all experimental results are known for really three-dimensional bubbles, a qualitative comparison is possible. As a main reference the book by [START_REF] Clift | Bubbles, Drops and Particles[END_REF] may be cited, where a general diagram of bubble shapes in dependence on the Eötvös and Reynolds numbers can be found (see Figure 7). The comparison enables us to conclude that our numerical bubble shapes are in a good agreement with the experimental predictions. Our results depicted in Figure 5 also compare well with the numerical results of [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF] (Figure 3 in [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF]), and the results of Figure 6 with the experimental results of [START_REF] Bhaga | Bubbles in viscous liquids: shapes, wakes and velocities[END_REF] (Figure 3e in [START_REF] Bhaga | Bubbles in viscous liquids: shapes, wakes and velocities[END_REF]; see also the numerical results of [START_REF] Chen | The development of a bubble rising in a viscous liquid[END_REF], Figures 15 and19).

The typical bubble shapes and the streamlines in the local reference frames of the bubbles are shown in Figure 8. It is apparent that all basic bubble shape regimes are successfully recovered with the physical parameters lying exactly within the experimentally predicted ranges (see Figure 7). It is interesting to notice that, while the secondary vortices develop inside of the bubble in any case, the vortices in the bubble wake are pronounced only when the bubble deformation is significant, i.e. for high Eötvös numbers (especially remarkable fact is the formation of the recirculating region behind the bubble in dimpled ellipsoidal cap and spherical cap regimes, see Figures 8c and8e). The results agree very well with the computations of [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF] and with the experiments of [START_REF] Bhaga | Bubbles in viscous liquids: shapes, wakes and velocities[END_REF].

Figure 9 illustrates the case of the so-called "wobbling bubble" regime. The wobbling typically appears with sufficiently high Reynolds numbers when the Eötvös number is, roughly, in the range between 1 and 100 (see Figure 7). In this shape-regime, the bubble initially has a nearly spherical cap shape; however, at later stage of the motion, a remarkable flattening of the bubble top can be observed (see Fig. 9). The bubble bottom undergoes permanent deformations resulting from the unstable and unsymmetric evolution of the bubble wake. In particular, the unsymmetric pairs of secondary vortices are clearly observed in the wake as the consequence of asynchronous separation of the boundary layer from different sides of the bubble surface. This flow pattern bears some resemblance to the von Karman vortex path typically formed behind a rigid body in a highly convective flow. With higher Eötvös numbers, one can observe even more pronounced wake instability (see Figure 10). While the wobbling of the bubble bottom remains apparent, the flattening of the bubble top is not seen, which makes the bubble shape resemble a spherical cap. The remarkable bubble wake, known in the literature as helical vortex path, is in a good qualitative agreement with experimental investigations of [START_REF] Tsuchiya | Near-wake structure of a single gas bubble in a two-dimensional liquid-solid fluidized bed: vortex shedding and wake size variation[END_REF].

The results of Figures 9 and 10 agree qualitatively well with the experimental observations cited in [START_REF] Clift | Bubbles, Drops and Particles[END_REF]; we are unaware of any other numerical simulations of a wobbling bubble dynamics.

Remark 4.1. The extensive numerical experiments have been conducted in the work of Smolianski et al. (2003a), where the computed terminal rise velocities as well as the positions of bubbles in dependence on time for diverse shape regimes have been shown to agree very well with the available experimental data.

The next group of numerical experiments is devoted to the convergence study. We take the dimensionless parameters as ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , Re = 35, Eo = 125 and compare the solutions corresponding to h = 1/20, 1/40, 1/80 at the time t = 3.0 s. Figure 11 illustrates the shape convergence, while Figure 12 shows the convergence of the pressure (the modified pressure, i.e. the pressure without its hydrostatic component, is depicted). We may conclude that with the present method the qualitatively correct picture may be obtained even on very coarse meshes with h ≈ 1/20. It is worth noting that the pressure discontinuity is captured very well; this fact has been already emphasized in the preceding section.

The following set of numerical tests is related to the mass conservation issue. We take again the physical parameters ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , Re = 35, Eo = 125, the 1 elements for velocity and the triangulation based on 20×40 rectangular grid. At the time t = 3.0s we compare the result obtained without reinitialization and level-set correction steps with the one computed with reinitialization-correction procedure. Figure 13 shows the relative change of bubble mass (i.e. of bubble area) as a function of time.

The necessity of the reinitialization-correction procedure is clearly confirmed by Figure 13. The maximal relative error in the bubble mass on the considered time interval [0; 3] is about 25% without reinitializationcorrection and 0.8% with the reinitialization-correction procedure.

Figure 14 illustrates the grid convergence for the relative change of bubble mass when the reinitialization and correction steps are performed. The maximal relative error in the bubble mass on the time interval [0; 3] is about 3h 2 for each of the three computations corresponding to h = 1/20, 1/40, 1/80. This confirms the theoretically predicted O(h 2 )-accuracy of mass conservation. The conservation of bubble mass is, indeed, very good: for h = 1/40 the maximal relative error is about 0.2%, and for h = 1/80 it is ≈ 0.05%.

Breaking bubble

In this test we use the same initial and boundary conditions as in the preceding section; the dimensionless parameters are ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , Re = 700, Eo = 500. For such Eötvös number the bubble should become "skirted" in accordance with the experimental predictions (see Figure 7), which can be also seen in Figure 15. Later, the skirt breaks off due to the action of the vortices in the bubble wake, and the remaining part of the bubble rapidly develops a spherical-cap shape. This is also in a full agreement with experimental observations (see [START_REF] Hnat | Spherical cap bubbles and skirt formation[END_REF]), but, in reality, it happens at much smaller Reynolds numbers (about 10÷100) than those used in our experiments. That is, probably, due to the essentially three-dimensional behavior of a real skirted bubble (see also [START_REF] Sussman | Axisymmetric free boundary problems[END_REF]). However, a qualitative picture of bubble break-up seems to be captured well even with the presented two-dimensional computations.

Our results compare well with the 2D numerical results of [START_REF] Sussman | A level set approach for computing solutions to incompressible twophase flow[END_REF], who employed the finitedifference level-set method, and with the 2D numerical predictions of [START_REF] Baker | The rise and distortion of a two-dimensional gas bubble in an inviscid liquid[END_REF], who advocated the boundary integral method for an inviscid gas bubble.

A typical time-dependence for the relative change of bubble mass is shown in Figure 16. The fluctuation in about 5% appears precisely at the moment of the bubble break-up and seems to be unavoidable, because of a limited grid-resolution inherent to any numerical experiment. However, the method quickly recovers the bubble mass, and the relative error of mass conservation becomes a small fraction of one percent.

Merger of two bubbles

Here we consider the rectangular domain (0; 1) × (0; 2) with two circular bubbles inside; the center of the first bubble is (0.5; 1.0) and its radius is equal to 0.25, the center of the second bubble is (0.5; 0.5) and the radius is 0.2. Thus, the bubbles have a common axis of symmetry, and the initial distance between them equals 1/5 of the radius of the largest bubble. We take zero velocity field at the initial moment and the same boundary conditions as in the preceding sections. The dynamics of the bubbles, to a large extent, depends on the initial distance between them and on the magnitude of the surface tension. If the surface tension is high enough (correspondingly, the Eötvös number is small), no merger happens, the bubbles develop nearly ellipsoidal shapes and rise separately (see, e.g., [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF]). Hence, in order to simulate a merger process, we take comparably small surface tension coefficient and, to get the merger earlier, a small initial distance between the bubbles. Our non-dimensional parameters for this test are Re = 35, Eo = 250, ρ 1 /ρ 2 = 10 2 , µ 1 /µ 2 = 10 (Reynolds and Eötvös numbers are based on the diameter of the large bubble).

Figure 17 illustrates the merger process. Since the small bubble is located very close to the large one, this lower bubble turns out to be in the wake of its upper "neighbour" and rises faster than that. In the process, two opposite signed vortices are created in the wake of the large bubble. This produces a lower pressure region behind the large bubble and generates flow streaming into the symmetry line of the flow. As a result, the front portion of the small bubble becomes narrower and sharper. At time t = 0.525, we see that the head of the small bubble almost catches up with the bottom of the large bubble. In the next moment t = 0.6, the two bubbles merge into a single bubble. At this time, the interface conjunction forms a cusp singularity that is rapidly smoothed out by viscosity and surface tension.

The results compare well with the computations of [START_REF] Chang | A level set formulation of Eulerian interface capturing methods for incompressible fluid flows[END_REF], [START_REF] Tornberg | Interface Tracking Methods with Application to Multiphase Flows[END_REF] and [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF]. The results compare favorably with the numerical predictions of [START_REF] Delnoij | Computational fluid dynamics (CFD) applied to dispersed gas-liquid two-phase flows[END_REF] who found a qualitative agreement with available experimental data.

The dependence on time for the relative change of total bubble mass is shown in Figure 18. The fluctuations in about 2% appear at the moment of the bubble merger. However, owing to the reinitialization-correction procedure, the method quickly recovers the total bubble mass and keeps the error of mass conservation within the predicted O(h 2 )-accuracy.

Rayleigh-Taylor instability

The Rayleigh-Taylor instability, associated with the acceleration of a heavy fluid into a light one under the action of a gravitational field, is generic to a wide range of physical phenomena, and many numerical simulations have been performed (see e.g. [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF][START_REF] Tryggvason | Computations of three-dimensional Rayleigh-Taylor instability[END_REF]). [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF] analyzed the problem by means of the linear theory, which is applicable in the early stages of instability development.

In our simulations we took the data as in [START_REF] Puckett | A high-order projection method for tracking fluid interfaces in variable density incompressible flows[END_REF] (the same data were used by [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF]). Namely, the computational domain is the rectangle (0; 1) × (0; 4), the viscosities of both fluids are equal to µ 1 = µ 2 = 3 • 10 -3 kg/(m • s), the fluid densities are ρ 2 = 1.2 kg/m 3 and ρ 1 = 0.17 kg/m 3 , and the surface tension coefficient σ is zero. The initial interface shape is given by the function y(x) = 2.0 + 0.05 cos(2πx) (see Figure 19), initial velocity field is zero, and the boundary conditions are no-slip on the top and the bottom walls and free-slip on the vertical walls.

Figure 20 shows the development of the instability. As the heavy fluid penetrates the light fluid, the interface begins to roll up along the sides of the spike giving the characteristic mushroom shape. This phenomenon, known as the Kelvin-Helmholtz instability, is due to the development of short wavelength perturbations along the fluid interface and parallel to the main flow. Our numerical predictions compare well with those of [START_REF] Puckett | A high-order projection method for tracking fluid interfaces in variable density incompressible flows[END_REF] and of [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF]. The maximum mass fluctuation is approximately 0.06%, which is better than in front-tracking based simulations of [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF] but slightly worse than in volume-of-fluid based computations of [START_REF] Puckett | A high-order projection method for tracking fluid interfaces in variable density incompressible flows[END_REF] (although, our grid was 1.5 times coarser).

Figure 21 illustrates the instability development in the presence of small amount of surface tension. The surface tension coefficient was chosen to be smaller than the critical value for which the flow is stabilized. According to [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF], for an initial perturbation whose wavelength is λ, this critical value is given by

σ c = λ 2 (ρ 2 -ρ 1 )g 4π 2 .
In our case λ = 1, hence σ c ≈ 2.6 • 10 -2 ; we took σ = 0.015. The regularizing effect of surface tension is clearly seen from Figure 21: the surface tension decreases the growth rate of the Rayleigh-Taylor instability, delays the development of the secondary, Kelvin-Helmholtz instability, and makes the interface smoother and more compact. This is in a good agreement with numerical observations of [START_REF] Lock | Local mesh adaptation technique for front tracking problems[END_REF].

Bifurcating jet

This test concerns the high-speed propagation of a horizontal viscous jet in a viscous medium. The jet evolution is considered in the absence of gravity, but in the presence of surface tension. The geometry consists of the rectangular domain (0; 3) × (0; 1) and the jet inflow boundary {(x 1 ; x 2 )| x 1 = 0, 0.45 ≤ x 2 ≤ 0.55}. The initial velocity field is zero; the inflow velocity is horizontal and its magnitude evolves in time as 5(1 -e -2t ).

The prescribed boundary conditions are the outflow (zero tangential and normal stress) at the right side, the inflow at the jet inflow part of the boundary, and the no-slip on the rest of the domain's boundary. We assume that both fluids have equal viscosity µ. Then, the problem is characterized by three nondimensional parameters: the density ratio ρ 2 /ρ 1 , the Reynolds number Re = ρ 2 U d/µ and the Weber number W e = ρ 1 U 2 d/σ, where index "1" corresponds here to the fluid around the jet, index "2" to the jet fluid, U is the magnitude of the inflow velocity and d is the inflow diameter. According to the values of these parameters, diverse regimes of jet flow can be observed (see [START_REF] Lin | Drop and spray formation from a liquid jet[END_REF]). Figure 22 illustrates the jet evolution with ρ 2 /ρ 1 = 10, Re = 5000 and W e = 250, which approximately corresponds to the so-called second windinduced regime. This is a high jet-velocity regime with relatively small surface tension coefficient. As the flow starts, the interface front expands very quickly and undergoes the interfacial stress fluctuations which cause the unstable growth of short-length waves. Those waves could be suppressed by the surface tension, but the interface is nearly flat in the vicinity of its frontal point, hence, the surface tension (that is proportional to the interface curvature) becomes almost negligible there; this results in the jet bifurcation. In our computations the jet break-up happens approximately at time t = 0.3 s at the distance 0.5 = 5d downstream, which is in excellent agreement with the predictions of [START_REF] Danaila | Direct numerical simulation of bifurcating jets[END_REF]. Generally, the jet unbroken length is a function of the dimensionless parameters of the problem. At later time, the jet branches start to roll up due to the Kelvin-Helmholtz instability, and multiple secondary bifurcations as well as drop formations can be observed (see Figure 22, bottom row). The qualitative behavior of the jet in our numerical simulation compares satisfactorily with experimental results, see [START_REF] Lin | Drop and spray formation from a liquid jet[END_REF].

Conclusions

A unified approach for the numerical simulation of unsteady two-phase flows with free moving interfaces has been presented. It essentially relies on the finite-element spatial approximation, the level-set interface representation and the operator-splitting temporal discretization. Such a combination of the numerical techniques results in a very simple and cost-effective algorithm that allows to simulate virtually all kinds of two-fluid interfacial flows (bubbles, jets, waves, films etc.). Despite its computational simplicity, the proposed approach demonstrates a good accuracy, which has been shown in diverse numerical tests (see also Smolianski et al. (2003a) and Smolianski et al. (2003b) where the present method has been used for a detailed study of bubble dynamics and numerous comparisons with experimental data have been done).

Unlike many existing numerical methods, the proposed approach is flexible enough to permit the simulations of diverse flow regimes ranging from the Stokes to highly convective flows. It allows for the modelling of flows with breaking and merging interfaces, while maintaining a sharp interface approximation throughout the computational process. The combination of the finite-element discretization and the level-set method makes it possible to model the interface without introducing any artificial numerical parameters like the interface thickness; the interface can be localized precisely, which enables us to accurately compute the surface tension force and the terms with discontinuous density/viscosity coefficient in the discrete momentum equation. We have also presented simple techniques allowing to obtain a 2nd-order accurate curvature approximation and to maintain 2nd-order accurate global mass balance. The numerical experiments show that the accuracy of the interface approximation and of the capillary force evaluation has a dominating role in the simulations of interfacial flows. This, possibly, explains a good performance of the presented approach, in spite of the fact that the method is formally only 1st-order accurate in time and in space. We have, in fact, compared the results with those obtained by a 2nd-order-in-space method (this can be readily achieved via the local mesh adjustment near the interface, see e.g. [START_REF] Lock | Local mesh adaptation technique for front tracking problems[END_REF]). The difference in accuracy for the velocity field appears to be very small with the realistic mesh sizes, and almost no difference can be seen in the interface shapes.

Finally, we may note that the method utilizes the simplest, piecewise-linear continuous approximation for all unknowns (velocity, pressure, interface) on a single uniform grid which is never changed in the process of computations and that the whole algorithm can be easily parallelized owing to the nature of operator-splitting approach. The work on the extension of the method to 3D case is currently under way, and the corresponding results will be reported in forthcoming papers.

Figure 1 :

 1 Figure 1: Sketch of a two-fluid flow configuration.

  characteristic time-scale should be such that all physical effects forming corresponding non-dimensional criterion counterbalance each other; thus, from F r = The surface tension phenomenon in viscous flow is characterized by the capillary number Ca = µU(ca) σ

Figure 3 :

 3 Figure 3: Spurious currents in a static bubble simulation.

  Ca 9.1 • 10 -3 9.0 • 10 -3 8.8 • 10 -3 9.1 • 10 -3 9.4 • 10 -3

Figure 4 :

 4 Figure 4: Pressure field corresponding to three different mesh sizes.

Figure 5 :

 5 Figure 5: Evolution of rising bubble with large surface tension (ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , Re = 35, Eo = 25); triangulation is based on 40 × 80 rectangular grid.

Figure 6 :

 6 Figure 6: Evolution of rising bubble with small surface tension (ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , Re = 35, Eo = 125); triangulation is based on 40 × 80 rectangular grid.

Figure 7 :

 7 Figure 7: Shape regimes for bubbles and drops. Figure reproduced from Clift et al. (1978) with permission of the publisher Academic Press.

Figure 8 :

 8 Figure 8: Different computed shapes of bubbles: (a) spherical with Re=1, Eo=0.6, (b) ellipsoidal with Re=20 Eo=1.2, (c) dimpled ellipsoidal cap with Re=35, Eo=125, (d) skirted with Re=55, Eo=875, (e) spherical cap with Re=94, Eo=115 and (f) wobbling with Re=1100, Eo=3.0; ρ 1 /ρ 2 = 10 3 ; µ 1 /µ 2 = 10 2 .

Figure 9 :

 9 Figure 9: The rise of a wobbling bubble (ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , Re = 1100, Eo = 3); triangulation is based on rectangular grid with h = 1/80.

Figure 10 :

 10 Figure 10: Helical vortex path behind the bubble (ρ 1 /ρ 2 = 10 2 , µ 1 /µ 2 = 10, Re = 900, Eo = 20.2); triangulation is based on rectangular grid with h = 1/80.

Figure 11 :

 11 Figure 11: Grid convergence test for the bubble shape (ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , Re = 35, Eo = 125); time t = 3.0.

Figure 12 :

 12 Figure 12: Grid convergence test for the modified pressure (ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , Re = 35, Eo = 125); time t = 3.0.

Figure 13 :

 13 Figure 13: Relative change of bubble mass vs. time (ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , Re = 35, Eo = 125; h = 1/20): "-• -" -without reinitialization-correction, "--" -with reinitialization-correction.

Figure 14 :

 14 Figure 14: Relative change of bubble mass vs. time (ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , Re = 35, Eo = 125; with reinitialization-correction): "--" -h = 1/20, "-• -" -h = 1/40, "--" -h = 1/80.

Figure 15 :

 15 Figure 15: Break-up of rising bubble; Re = 700, Eo = 500, ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , triangulation is based on rectangular grid with h = 1/80, 1 elements for velocity, ∆t = h/2.

Figure 16 :

 16 Figure 16: Breaking bubble: relative change of bubble mass vs. time; Re = 700, Eo = 500, ρ 1 /ρ 2 = 10 3 , µ 1 /µ 2 = 10 2 , h = 1/80.

Figure 17 :

 17 Figure 17: Merger of two rising bubbles; Re = 35, Eo = 250, ρ 1 /ρ 2 = 10 2 , µ 1 /µ 2 = 10, triangulation is based on rectangular grid with h = 1/40, 1 elements for velocity, ∆t = h/2.

Figure 18 :

 18 Figure 18: Merger of two bubbles: relative change of total bubble mass vs. time; Re = 35, Eo = 250, ρ 1 /ρ 2 = 10 2 , µ 1 /µ 2 = 10, h = 1/40.

Figure 19 :

 19 Figure 19: Initial configuration for Rayleigh-Taylor instability problem; the interface is given as y(x) = 2.0 + 0.05 cos(2πx).

Figure 20 :

 20 Figure 20: Rayleigh-Taylor instability with σ = 0; triangulation is based on rectangular grid with h = 1/40, 1 elements for velocity, ∆t = 2h.

Figure 21 :Figure 22 :

 2122 Figure 21: Rayleigh-Taylor instability with σ = 0.015; triangulation is based on rectangular grid with h = 1/40, 1 elements for velocity, ∆t = h.

Acknowledgements

The author would like to thank prof. Heikki Haario and Pasi Luukka from Technical University of Lappeenranta for the fruitful collaboration. This work was supported by the grant 70139/98 of Tekes, the National Technology Agency of Finland.