Advances in the simulation of multi-fluid flows with the Particle Finite Element Method. Application to bubble dynamics

Abstract

In this work we extend the Particle Finite Element Method (PFEM) to multi-fluid flow problems with the aim of exploiting the fact that Lagrangian methods are specially well suited for tracking interfaces. We develop a numerical scheme able to deal with large jumps in the physical properties, included surface tension, and able to accurately represent all types of discontinuities in the flow variables. The scheme is based on decoupling the velocity and pressure variables through a pressure segregation method which takes into account the interface conditions. The interface is defined to be aligned with the moving mesh, so that it remains sharp along time, and pressure degrees of freedom are duplicated at the interface nodes to represent the discontinuity of this variable due to surface tension and variable viscosity. Furthermore, the mesh is refined in the vicinity of the interface to improve the accuracy and the efficiency of the computations.

We apply the resulting scheme to the benchmark problem of a two-dimensional bubble rising in a liquid column presented in [1], and propose two breakup and coalescence problems to assess the ability of a multi-fluid code to model topology changes.

Keywords: Particle Finite Element Method (PFEM); Lagrangian simulation; multi-fluid flows; pressure segregation; surface tension; bubble dynamics.

1 INTRODUCTION

The simultaneous presence of multiple fluids with different properties in external or internal flows is found in daily life, environmental problems, and numerous industrial processes. Examples are fluid-fuel interaction in enhanced oil recovery, blending of polymers, emulsions in food manufacturing, rain droplet formation in clouds, fuel injection in engines, and bubble column reactors, to name only a few. Although multi-fluid flows occur frequently in nature and engineering practice, they still pose a major research challenge from both theoretical and computational points of view.

In the case of immiscible fluids, the dynamics of the interface between fluids plays a dominant role. The success of the simulation of such flows depends on the ability of the numerical method to model accurately the interface and the phenomena taking place on it, such as the surface tension. The origin of the surface tension force lies in the different intermolecular attractive forces that act in the two fluids, and the result is an interfacial energy per area that acts to resist the creation of new interface, so that the interface behaves like a stretched membrane trying to minimize its area.

The main difference between a single-fluid flow (homogeneous) and a multi-fluid (heterogeneous) one is the presence of internal interfaces. In addition to the well-known difficulties in the simulation of homogeneous flows, namely the coupling of pressure and velocity through the incompressibility constraint, the need of the discretization spaces to satisfy the inf-sup condition, and the non-linearity of the governing equations, numerical methods for heterogeneous flows face the following challenges:

1. Accurate definition of the interface position. The interface separating the fluids needs to be tracked accurately without introducing excessive numerical smoothing.
2. Modeling the jumps in the fluid properties across the interface. Large jumps of fluid density and viscosity across the interface need to be properly taken into account in order to satisfy the momentum balance at the vicinity of the interface.
3. Modeling the discontinuities of the flow variables across the interface. Velocity and pressure may be discontinuous across the interface under certain conditions.
4. Modeling the surface tension. Since surface tension plays a very important role in the immiscible interface dynamics, this force needs to be accurately evaluated and incorporated into the model.

This paper extends the work of the authors in the study of multi-fluid flows with the Particle Finite Element Method [2]. The new contributions include the scheme for describing the interface, the computation of surface tension and the enhanced pressure segregation approach. We have found that these developments are essential for accurately modeling of bubble dynamics. The paper is organized as follows: In Section 2 we present the governing equations together with the interface conditions and the possible discontinuities of the flow variables. In Section 3 we briefly review the interface descriptions most used in the literature, to focus later in the Particle Finite Element Method (Section 4) and the numerical scheme we have developed (Sections 5 and 6). Finally, in Section 7 we apply the PFEM to the solution of the two-dimensional bubble rise, breakup and coalescence problems.

2 GOVERNING EQUATIONS

Figure 1: Two-fluid flow configuration.

Let ${\textstyle \Omega \subset {\boldsymbol {R}}^{d}}$, ${\textstyle d\in \{2,3\}}$, be a bounded domain containing two different fluids (see Figure 1). We denote time by ${\textstyle t}$, the Cartesian spatial coordinates by ${\textstyle {\boldsymbol {x}}=\{x_{i}\}_{i=1}^{d}}$, and the vectorial operator of spatial derivatives by ${\textstyle {\boldsymbol {\nabla }}=\{\partial _{x_{i}}\}_{i=1}^{d}}$. The evolution of the velocity ${\textstyle u=u({\boldsymbol {x}},t)}$ and the pressure ${\textstyle p=p({\boldsymbol {x}},t)}$ is governed by the Navier-Stokes equations:

$\rho \displaystyle {\frac {{d}u}{{d}t}}={\boldsymbol {\nabla }}\cdot {\boldsymbol {\sigma }}+\rho {\boldsymbol {g}}\quad {\mbox{ in}}\;\Omega \times (0,T)$ (1.a) $\displaystyle {\frac {{d}\rho }{{d}t}}+\rho {\boldsymbol {\nabla }}\cdot {u}=0\quad {\mbox{ in}}\;\Omega \times (0,T)$ (1.b)

where ${\textstyle \rho }$ is the density, ${\textstyle {\boldsymbol {\sigma }}}$ the Cauchy stress tensor, ${\textstyle {\boldsymbol {g}}}$ the vector of gravity acceleration, and ${\textstyle \displaystyle {\frac {{d}\phi }{{d}t}}}$ represents the total or material derivative of a function ${\textstyle \phi }$. The constitutive equation for a Newtonian and isotropic fluid takes the form

${\boldsymbol {\sigma }}=-p{\boldsymbol {I}}+2\mu ({\boldsymbol {D}}-{\frac {1}{3}}\varepsilon _{V}{\boldsymbol {I}})$

(2)

with ${\textstyle {\boldsymbol {I}}}$ the identity tensor, ${\textstyle \mu }$ the dynamic viscosity, ${\textstyle {\boldsymbol {D}}={\frac {1}{2}}({\boldsymbol {\nabla }}u+{\boldsymbol {\nabla }}^{T}u)}$ the strain rate tensor, and ${\textstyle \varepsilon _{V}={\boldsymbol {\nabla }}\cdot {u}}$ the volumetric strain rate.

Let ${\textstyle \Gamma _{int}(t)}$ be the interface that cuts the domain ${\textstyle \Omega }$ in two open subdomains, ${\textstyle \Omega {^{+}}(t)}$ and ${\textstyle \Omega {^{-}}(t)}$, which satisfy: ${\textstyle \Omega {^{+}}\cap \Omega {^{-}}=\emptyset }$, ${\textstyle \Omega ={\bar {\Omega }}{^{+}}\cup {\bar {\Omega }}{^{-}}}$, and ${\textstyle \Gamma _{int}={\bar {\Omega }}{^{+}}\cap {\bar {\Omega }}{^{-}}=\partial \Omega {^{+}}\cap \partial \Omega {^{-}}}$. In each subdomain, the physical properties are defined as:

$\rho =\rho ({\boldsymbol {x}},t)=\left\{{\begin{array}{ll}\rho {^{+}}&{\mbox{ if }}{\boldsymbol {x}}\in \Omega {^{+}}(t)\\\rho {^{-}}&{\mbox{ if }}{\boldsymbol {x}}\in \Omega {^{-}}(t)\end{array}}\right.,\qquad \mu =\mu ({\boldsymbol {x}},t)=\left\{{\begin{array}{ll}\mu {^{+}}&{\mbox{ if }}{\boldsymbol {x}}\in \Omega {^{+}}(t)\\\mu {^{-}}&{\mbox{ if }}{\boldsymbol {x}}\in \Omega {^{-}}(t)\end{array}}\right.$

(3)

If density and viscosity are assumed to remain constant in each fluid (i.e. fluids are incompressible, immiscible, and isothermal), we have that ${\textstyle \displaystyle {\frac {{d}\rho }{{d}t}}=0}$ and ${\textstyle \displaystyle {\frac {{d}\mu }{{d}t}}=0}$. Consequently, we have on the one side that ${\textstyle \varepsilon _{V}={\boldsymbol {\nabla }}\cdot {u}=0}$, this is the mass conservation equation for incompressible flows; and on the other side, that ${\textstyle \displaystyle {\frac {d}{{d}t}}\,\Gamma _{int}=0}$. This latter consequence means that interfaces are material surfaces, which move with the fluid velocity ${\textstyle u}$, and therefore, they are naturally tracked in Lagrangian formulations.

Boundary and interface conditions

In order for the Navier-Stokes problem (1) to be well-posed, suitable boundary conditions need to be specified. On the external boundary ${\textstyle \partial \Omega =\Gamma _{D}\cup \Gamma _{N}}$, such that ${\textstyle \Gamma _{D}\cap \Gamma _{N}=\emptyset }$, we consider the following:

$u={\bar {u}}\quad {\mbox{ on}}\;\Gamma _{D}$ (4) ${\boldsymbol {\sigma }}\cdot {\boldsymbol {n}}={\bar {\boldsymbol {\sigma }}}_{n}\quad {\mbox{ on}}\;\Gamma _{N}$ (5)

${\textstyle {\bar {u}}}$ is the prescribed velocity, ${\textstyle {\boldsymbol {n}}}$ the outer unit normal to ${\textstyle \Gamma _{N}}$, and ${\textstyle {\bar {\boldsymbol {\sigma }}}_{n}}$ the prescribed traction vector. A Neumann boundary ${\textstyle \Gamma _{N}}$ with ${\textstyle {\bar {\boldsymbol {\sigma }}}_{n}={\boldsymbol {0}}}$ is called free surface.

On the internal interfaces ${\textstyle \Gamma _{int}}$, the coupling conditions are [3]:

${\mbox{⟦}}u{\mbox{⟧}}={\boldsymbol {0}}\qquad {\mbox{ on}}\;\Gamma _{int}$ (6) ${\mbox{⟦}}{\boldsymbol {\sigma }}{\mbox{⟧}}\cdot {\boldsymbol {n}}=\gamma \kappa {\boldsymbol {n}}\qquad {\mbox{ on}}\;\Gamma _{int}$ (7)

with ${\textstyle {\boldsymbol {n}}}$ now the unit normal to ${\textstyle \Gamma _{int}}$, ${\textstyle \gamma }$ the surface tension coefficient, ${\textstyle \kappa }$ the interface curvature, and ${\textstyle {\mbox{⟦}}\phi {\mbox{⟧}}=\phi {^{+}}-\phi {^{-}}}$ represents the jump of a quantity ${\textstyle \phi }$ across the interface.

Equation (6) expresses the continuity of all velocity components. The normal component has to be continuous when there is no mass flow through the interface, and the tangential components have to be continuous when both fluids are viscous (${\textstyle \mu ^{+},\mu ^{-}>0}$), similar to a no-slip condition.

Equation (7) expresses that the jump in the normal stresses is balanced with the surface tension force. This force is proportional to the interface curvature and points to the center of the osculating circle that approximates ${\textstyle \Gamma _{int}}$. The surface tension coefficient ${\textstyle \gamma }$ is assumed constant and its value depends on the two fluids at the interface.

2.1 Interface discontinuities

Discontinuities at the interface can be of two types:

• ${\cal {C}}^{0}$ discontinuity, when the flow variable has a kink (i.e. the gradient has a jump), and
• ${\cal {C}}^{-1}$ discontinuity, when the flow variable itself has a jump.

Differences in density at the interface cause a kink in the hydrostatic pressure profile, leading to a jump in the pressure gradient, and then to a ${\textstyle {\cal {C}}^{0}}$ discontinuity in the pressure field (Fig. 2a).

Differences in viscosity lead to discontinuous components of the strain rate tensor ${\textstyle {\boldsymbol {D}}}$, and therefore to a ${\textstyle {\cal {C}}^{0}}$ discontinuity of the velocity field at the interface (Fig. 2b):

${\boldsymbol {t}}\cdot {\mbox{⟦}}{\boldsymbol {\sigma }}{\mbox{⟧}}\cdot {\boldsymbol {n}}=0\quad \Longrightarrow \quad \mu ^{+}\left(\displaystyle {\frac {\partial u_{t}}{\partial n}}+\displaystyle {\frac {\partial u_{n}}{\partial s}}\right)^{+}-\mu ^{-}\left(\displaystyle {\frac {\partial u_{t}}{\partial n}}+\displaystyle {\frac {\partial u_{n}}{\partial s}}\right)^{-}=0$

(8)

with ${\textstyle \partial _{s}={\boldsymbol {t}}\cdot {\boldsymbol {\nabla }}}$ the tangential derivative.

Both differences in viscosity and the presence of surface tension cause a ${\textstyle {\cal {C}}^{-1}}$ discontinuity in the pressure field (Figs. 2b and 2c), as shown in [4]:

${\boldsymbol {n}}\cdot {\mbox{⟦}}{\boldsymbol {\sigma }}{\mbox{⟧}}\cdot {\boldsymbol {n}}=\gamma \kappa \quad \Longrightarrow \quad p^{+}-p^{-}=2(\mu ^{+}-\mu ^{-})\displaystyle {\frac {\partial u_{n}}{\partial n}}-\gamma \kappa$

(9)

Notice that even in the case of ${\textstyle \gamma =0}$, pressure is discontinuous when ${\textstyle \mu ^{+}\neq \mu ^{-}}$.

(a)	(b)
(c)
Figure 2: Flow discontinuities for: (a) density jump, (b) viscosity jump, and (c) surface tension.

3 INTERFACE DESCRIPTION

A major challenge in the simulation of interfaces between different fluids is the accurate description of the interface evolution. The location of the interface is in general unknown and coupled to the local flow field which transports the interface. It is essential that the interface remains sharp along time and is able to fold, break and merge. In the past decades a number of techniques have been developed to model interfaces in multi-fluid flow problems, each technique with its own particular advantages and disadvantages. Comprehensive reviews can be found in e.g. [5,6,7,8].

(a)	(b)
Figure 3: (a) Moving mesh adapted to the interface, and (b) fixed mesh, where interface moves through the elements.

The main classification of interface descriptions is regarding the reference frame adopted (see Figure 3). In the moving mesh methods, the mesh is deformable and adapted to the interface, which is explicitly tracked along the trajectories of the fluid particles. Examples are methods based on the Arbitrary Lagrangian-Eulerian (ALE) formulation [9,10], the deformable-spatial-domain/stabilized space-time deformation (DSD/SST) method [11,12], or the fully Lagrangian formulation such as in [13,14] and the Particle Finite Element Method [15,16,17,18].

On the other hand, fixed mesh methods use a separate procedure to describe the position of the interface. They can be further grouped in front-tracking methods, which use massless marker points to follow the fluid interface while the Navier-Stokes equations are solved on a fixed mesh [19,20], and front-capturing methods, which introduce a new variable ${\textstyle \psi }$ in the model to describe the presence or not of a fluid in a position of the domain. The most extended front-capturing methods are the Volume-Of-Fluid, originally developed by Hirt [21], and the Level Set method by Osher et al. [22].

4 PARTICLE FINITE ELEMENT METHOD

The Particle Finite Element Method (PFEM) is a Lagrangian numerical technique for modeling and analysis of complex multidisciplinary problems in fluid and solid mechanics involving thermal effects, interfacial and free-surface flows, and fluid-structure interaction, among others.

PFEM is a particle method in the sense that the domain is defined by a collection of particles that move in a Lagrangian manner according to the calculated velocity field, transporting their momentum and physical properties (e.g. density, viscosity). The interacting forces between particles are evaluated with the help of a mesh. Mesh nodes coincide with the particles, so that when the particles move so does the mesh. On this moving mesh, the governing equations are discretized using the standard finite element method. The possible large distortion of the mesh is avoided through an efficient Delaunay triangulation remeshing [23] of the computational domain. Due to the fact that all the hydrodynamical information is stored in the nodes, remeshing does not introduce numerical diffusion. This gives the method excellent capabilities for modeling large displacement and large deformation problems. A more detailed explanation of the method can be found in e.g. [16,17,24].

Since in PFEM the interface is described by mesh nodes and element edges (Fig. 3a), it is a well-defined curve, with the information regarding its location and curvature readily available. The interface nodes carry the jump of properties (e.g. density, viscosity), maintaining the interface sharp without diffusion along time. Furthermore, it is straightforward to impose the boundary conditions on the interface and to treat any number of fluids.

Regarding the modeling of the jumps in the fluid properties across the interface, while in fixed mesh methods typically the interface is considered to have a finite thickness and the fluid properties change smoothly and continuously from the value on the one side of the interface to the value on the other side, PFEM treats the interface in a sharp manner, so that it is clear which property value is valid at each point.

Figure 4: Pressure profiles when using continuous and discontinuous representations.

Regarding the modeling of the discontinuities in the flow variables across the interface (see Section 2.1), in fixed mesh methods where the physical properties have been smoothed, functions are continuous across the interface and thus not appropriate for the approximation of discontinuous variables. When the physical properties are modeled sharp, the elements cut by the interface require a special treatment in order to be able to represent the discontinuities. Gravity dominated flows will require “enrichment” of the pressure approximation, and viscosity dominated flows will require “enrichment” of the velocity approximation. On the contrary, ${\textstyle {\cal {C}}^{0}}$ discontinuities need no special attention when the interface is aligned with the mesh, as the kinks in the solution are automatically represented. Only ${\textstyle {\cal {C}}^{-1}}$ discontinuities need some attention in PFEM. In particular, the pressure field has been made double-valued at the interface, i.e. pressure degrees of freedom have been duplicated (${\textstyle p^{+}}$, ${\textstyle p^{-}}$) in the interface nodes [4]. The pressure discontinuity is thus optimally approximated. Figure 4 shows that the use of continuous pressure representations may introduce errors in the incompressibility condition.

5 PRESSURE SEGREGATION FOR THE MULTI-FLUID EQUATIONS

The Lagrangian approach simplifies the equations by separating the problem into a geometrical part (tracking the motion of the nodes)

$\displaystyle {\frac {{d}{\boldsymbol {x}}}{{d}t}}=u({\boldsymbol {x}},t)$

(10)

and a physical part (calculating how the flow variables change in time at each node) described by the following Stokes equations:

$\rho ({\boldsymbol {x}})\displaystyle {\frac {{d}u}{{d}t}}-{\boldsymbol {\nabla }}\cdot {2}\mu ({\boldsymbol {x}}){\boldsymbol {D}}(u)+{\boldsymbol {\nabla }}p=\rho ({\boldsymbol {x}}){\boldsymbol {g}}\quad {\mbox{in}}~\Omega$ (11.a) ${\boldsymbol {\nabla }}\cdot {u}=0\quad {\mbox{ in}}~\Omega$ (11.b) $u={\bar {u}}\quad {\mbox{ on}}~\Gamma _{D}$ (11.c) ${\boldsymbol {\sigma }}\cdot {\boldsymbol {n}}={\bar {\boldsymbol {\sigma }}}_{n}\quad {\mbox{ on}}\;\Gamma _{N}$ (11.d) ${\mbox{⟦}}u{\mbox{⟧}}=0\quad {\mbox{ and}}\quad {\mbox{⟦}}{\boldsymbol {\sigma }}{\mbox{⟧}}\cdot {\boldsymbol {n}}=\gamma \kappa {\boldsymbol {n}}\quad {\mbox{ on}}~\Gamma _{int}$ (11.e) $u({\boldsymbol {x}},0)=u^{0}({\boldsymbol {x}})\quad {\mbox{ in}}~\Omega$ (11.f)

The latter can be further split into two parts: one related to viscosity effects, and the other related to the incompressibility. According to the incremental pressure segregation scheme [25], and after implicit backward Euler time discretization, we propose the following splitting of equations (11) [26]:

1. Find an intermediate velocity ${\textstyle {\tilde {\boldsymbol {u}}}}$ solution of

   $\rho ({\boldsymbol {x}}){\frac {{\tilde {\boldsymbol {u}}}-{\boldsymbol {u}}^{n}}{\Delta t}}-{\boldsymbol {\nabla }}\cdot {2}\mu ({\boldsymbol {x}}){\boldsymbol {D}}({\tilde {\boldsymbol {u}}})+{\boldsymbol {\nabla }}p^{n}=\rho ({\boldsymbol {x}}){\boldsymbol {g}}$ (12.a) ${\tilde {\boldsymbol {u}}}={\bar {u}}\quad {\mbox{ on}}~\Gamma _{D}$ (12.b) $(-p^{n}{\boldsymbol {I}}+2\mu ({\boldsymbol {x}}){\boldsymbol {D}}({\tilde {\boldsymbol {u}}}))\cdot {\boldsymbol {n}}={\bar {\boldsymbol {\sigma }}}_{n}\quad {\mbox{ on}}~\Gamma _{N}$ (12.c) ${\mbox{⟦}}{\tilde {\boldsymbol {u}}}{\mbox{⟧}}=0\quad {\mbox{ and}}\quad {\mbox{⟦}}-p^{n}{\boldsymbol {I}}+2\mu ({\boldsymbol {x}}){\boldsymbol {D}}({\tilde {\boldsymbol {u}}}){\mbox{⟧}}\cdot {\boldsymbol {n}}=\gamma \kappa {\boldsymbol {n}}\quad {\mbox{ on}}~\Gamma _{int}</li>$ (12.d)

2. Determine ${\textstyle p^{n+1}}$ as solution of

   ${\boldsymbol {\nabla }}\cdot {\frac {1}{\rho ({\boldsymbol {x}})}}{\boldsymbol {\nabla }}(p^{n+1}-p^{n})={\frac {1}{\Delta t}}{\boldsymbol {\nabla }}\cdot {\tilde {\boldsymbol {u}}}\quad {in}~\Omega ,$ (13.a) ${\mbox{ and}}\quad {\boldsymbol {n}}\cdot {\boldsymbol {\nabla }}(p^{n+1}-p^{n})=0\quad {on}~\Gamma _{D}$ (13.b) $(p^{n+1}-p^{n}){\boldsymbol {n}}={\boldsymbol {0}}\quad {\mbox{ on}}~\Gamma _{N}$ (13.c) ${\mbox{⟦}}p^{n+1}-p^{n}{\mbox{⟧}}=0\quad {\mbox{ on}}~\Gamma _{int}</li>$ (13.d)

3. Finally, the ${\textstyle {\boldsymbol {u}}^{n+1}}$ velocity is obtained by

   ${\boldsymbol {u}}^{n+1}={\tilde {\boldsymbol {u}}}-{\frac {\Delta t}{\rho ({\boldsymbol {x}})}}{\boldsymbol {\nabla }}(p^{n+1}-p^{n})</li>$

   (14)

The interfacial stress jump condition has been consistently split between steps 1 and 2. Although the continuity of the pressure difference ${\textstyle p^{n+1}-p^{n}}$ across the interface expressed in Eq. (13.d) seems to be in contradiction with the fact that pressure is discontinuous at the interface, the scheme is able to capture this discontinuity without need of enforcing it explicitly [26].

After discretization in space (Galerkin weighted residual method), and assuming ${\textstyle \Gamma _{N}}$ to be a free surface (i.e. ${\textstyle {\bar {\boldsymbol {\sigma }}}_{n}={\boldsymbol {0}}}$), the split discrete equations read:

$1.\qquad {\frac {{\boldsymbol {M}}_{\rho }}{\Delta t}}({\tilde {\boldsymbol {U}}}-{\boldsymbol {U}}^{n})+{\boldsymbol {\boldsymbol {K}}}_{\mu }{\tilde {\boldsymbol {U}}}-{\boldsymbol {B}}{\boldsymbol {P}}^{n}={\boldsymbol {F}}^{n}$ (15) $2.\qquad \Delta t\,{\boldsymbol {\boldsymbol {L}}}_{(1/\rho )}{\boldsymbol {P}}^{n+1}=-{\boldsymbol {B}}^{T}{\tilde {\boldsymbol {U}}}+\Delta t\,{\boldsymbol {\boldsymbol {L}}}_{(1/\rho )}{\boldsymbol {P}}^{n}$ (16) $3.\qquad {\frac {{\boldsymbol {M}}_{\rho }}{\Delta t}}({\boldsymbol {U}}^{n+1}-{\tilde {\boldsymbol {U}}})-{\boldsymbol {B}}({\boldsymbol {P}}^{n+1}-{\boldsymbol {P}}^{n})={\boldsymbol {0}}$ (17)

where ${\textstyle {\boldsymbol {U}}}$, ${\textstyle {\boldsymbol {P}}}$ are the vectors of nodal velocities and pressure, ${\textstyle {\tilde {\boldsymbol {U}}}}$ is the intermediate velocity introduced by the fractional step, and ${\textstyle \Delta t}$ the time step. The matrices ${\textstyle {\boldsymbol {M}}_{\rho }}$ density weighted mass matrix, ${\textstyle {\boldsymbol {B}}}$ gradient matrix, ${\textstyle -{\boldsymbol {B}}^{T}}$ divergence matrix, ${\textstyle {\boldsymbol {\boldsymbol {K}}}_{\mu }}$ viscosity weighted stiffness matrix and the external force vector are defined as

${\boldsymbol {M}}_{\rho }^{ab}=\int _{\Omega }{\boldsymbol {N}}^{a}\rho {\boldsymbol {N}}^{b}\;d\Omega$ (18) ${\boldsymbol {\boldsymbol {K}}}_{\mu }^{ab}=\int _{\Omega }{\boldsymbol {\nabla }}{\boldsymbol {N}}^{a}\mu ({\boldsymbol {\nabla }}{\boldsymbol {N}}^{b}+{\boldsymbol {\nabla }}^{T}{\boldsymbol {N}}^{b})\;d\Omega$ (19) ${\boldsymbol {B}}^{ab}=\int _{\Omega }({\boldsymbol {\nabla }}\cdot {\boldsymbol {N}}^{a})\;N^{b}\;d\Omega$ (20) ${\boldsymbol {\boldsymbol {L}}}_{(1/\rho )}^{ab}=\int _{\Omega }{\boldsymbol {\nabla }}N^{a}{\frac {1}{\rho }}{\boldsymbol {\nabla }}N^{b}\;d\Omega$ (21) ${\boldsymbol {F}}^{a}=\int _{\Omega }{\boldsymbol {N}}^{a}\rho {\boldsymbol {g}}\,d\Omega +\int _{\Gamma _{int}}{\boldsymbol {N}}^{a}\gamma \kappa {\boldsymbol {n}}\,d\Gamma$ (22)

where ${\textstyle {\boldsymbol {N}}}$ and ${\textstyle N}$ are the standard linear FE shape functions for velocity and pressure, and superscripts ${\textstyle a,b}$ refer to node indices.

The accuracy in the treatment of the discontinuous density and viscosity will be determined by how well the numerical method is able to evaluate the integrals where these discontinuities are included. Unless the interface coincides with edges of elements, there is no way for standard element shape functions to capture the discontinuity of properties inside an element. This implies that one has either to increase the number of Gauss points (e.g. see discussion about numerical integration of discontinuous and singular functions in [27]) or enrich the shape function space, as e.g. in [28,29] and in the eXtended Finite Element Method (XFEM) [30,31]. The nodal interface description for immiscible fluids in PFEM allows to evaluate exactly these integrals.

Conditions (12.c) and (12.d) are naturally included in ${\textstyle {\boldsymbol {F}}}$, while pressure conditions (13.c) on ${\textstyle \Gamma _{N}}$ and (13.d) on ${\textstyle \Gamma _{int}}$ need to be weakly imposed in order to overcome the singularity of ${\textstyle {\boldsymbol {\boldsymbol {L}}}_{(1/\rho )}}$:

$\int _{\Gamma _{N}}q\,(p^{n+1}-p^{n})\,d\Gamma =0\qquad {\mbox{and}}\qquad \int _{\Gamma _{int}}q\,{\mbox{⟦}}p^{n+1}-p^{n}{\mbox{⟧}}\,d\Gamma =0$

Both integrals are discretized in space as ${\textstyle {\boldsymbol {M}}_{c}({\boldsymbol {P}}^{n+1}-{\boldsymbol {P}}^{n})}$, with ${\textstyle {\boldsymbol {M}}_{c}^{ab}=\int _{\Gamma }N^{a}N^{b}\,d\Gamma }$ the pressure mass matrix on the contours ${\textstyle \Gamma _{N}}$ and ${\textstyle \Gamma _{int}}$, and incorporated into the pressure Poisson equation in the following way [32]:

$\left(\Delta t\,{\boldsymbol {\boldsymbol {L}}}_{(1/\rho )}+\lambda \,{\boldsymbol {M}}_{c}\right){\boldsymbol {P}}^{n+1}=-{\boldsymbol {B}}^{T}{\tilde {\boldsymbol {U}}}+\left(\Delta t\,{\boldsymbol {\boldsymbol {L}}}_{(1/\rho )}+\lambda \,{\boldsymbol {M}}_{c}\right){\boldsymbol {P}}^{n}$

(23)

The penalty parameter ${\textstyle \lambda }$ is defined as ${\textstyle \lambda =\alpha \;{\frac {\Delta t}{h}}}$. It has to be sufficiently large so that the system matrix becomes invertible. ${\textstyle \alpha }$ is a scalar factor ${\textstyle \alpha \sim {\cal {O}}(1)}$ such that for ${\textstyle \alpha \rightarrow 0}$, the pressure equation is satisfied but the discrete system may continue singular, and ${\textstyle \alpha \rightarrow \infty }$ is equivalent to apply the pressure condition strongly and then the equation may not be satisfied.

Moreover, stabilization is needed in incompressible flows when interpolation spaces for velocity and pressure do not satisfy the inf-sup condition. Many stabilization procedures have been proposed in the literature, such as the Streamline-Upwind/Petrov-Galerkin, Galerkin Least-Squares, Finite Calculus or Orthogonal Sub-Scale methods. Those that include a projection of the pressure gradient need to be modified when density changes at the interface to take into account the ${\textstyle {\cal {C}}^{-1}}$ discontinuity of the hydrostatic pressure gradient. Details are explained in [2] and [33].

6 SURFACE TENSION

The surface tension force at the interface, ${\textstyle {\boldsymbol {f}}_{st}=\int _{\Gamma _{int}}\gamma \kappa {\boldsymbol {n}}\cdot {\boldsymbol {w}}\,d\Gamma }$, is naturally incorporated in the weak form of the finite element method. If it is discretized explicitly, i.e. the surface tension force is evaluated on the interface at the previous time step, the stability of the scheme imposes the following restriction on the time step size [34]:

$\Delta t_{st}<{\sqrt {\frac {\langle \rho \rangle \,h^{3}}{2\pi \gamma }}}$

(24)

where ${\textstyle \langle \rho \rangle ={\frac {1}{2}}(\rho _{1}+\rho _{2})}$, ${\textstyle h}$ is the mesh size and ${\textstyle \Delta t_{st}}$ the capillary time step. With this restriction the propagation of capillary waves is resolved and their unstable amplification avoided. Unfortunately, Eq. (24) can be rather limiting for fine meshes and large surface tension coefficients. An implicit [35] or semi-implicit [36] treatment of the surface tension would circumvent this constrain.

The interface representation by nodes and element edges in PFEM allows for an easy and accurate incorporation of the surface tension, avoiding the need of regularization techniques such as the Continuum Surface Force [34].

6.1 Curvature calculation

There are several ways to calculate the curvature ${\textstyle \kappa }$ from the information of the interface location. The one we have followed in this work is based on the osculating circle of a curve, which is defined as the circle that approaches the curve most tightly among all tangent circles at a given point. From the radius of the osculating circle, the quantity ${\textstyle \kappa {\boldsymbol {n}}}$ required for ${\textstyle {\boldsymbol {f}}_{st}}$ is calculated as:

${\boldsymbol {n}}={\frac {\boldsymbol {R}}{|{\boldsymbol {R}}|}},\qquad \kappa {\boldsymbol {n}}={\frac {\boldsymbol {R}}{|{\boldsymbol {R}}|^{2}}}$

(25)

Figure 5: Calculation of the osculating circle at node ${\boldsymbol {x}}$.

6.2 Static bubble

We verify our surface tension algorithm with the static bubble test case [37,26]. It consists in a circular fluid bubble into another viscous fluid at rest (see Figure 6), where gravitational or other external forces are neglected. According to the Laplace-Young law, the pressure jump will be ${\textstyle p_{2}-p_{1}=\gamma \kappa =\gamma /R}$, where ${\textstyle p_{2}}$ is the bubble internal pressure, ${\textstyle p_{1}}$ the outer pressure and ${\textstyle R}$ the bubble radius. Even with non-zero surface tension, the circular shape of the bubble should be preserved and the fluids should remain at rest no matter how long we integrate the equations in time.

Figure 6: Static bubble initial configuration.

We have simulated the equilibrium state of a circular bubble of a radius ${\textstyle R=0.25}$ (constant curvature ${\textstyle \kappa =1/R=4}$) in a static fluid with ${\textstyle \rho _{1}=\rho _{2}=1}$, ${\textstyle \mu _{1}=\mu _{2}=1}$, ${\textstyle g=0}$, and ${\textstyle \gamma =1}$. The simulations have been run 100 time steps with ${\textstyle \Delta t}$ fixed to 0.01. The bubble should remain exactly stationary and the velocity of the fluid should be exactly zero. But if the pressure is approximated continuously, the pressure fluctuations near the interface generate spurious velocity currents (see Figure 7) that may deform the bubble shape, produce a significant mass loss [38] and spoil the solution.

Figure 7: Spurious velocities in static bubble for continuous pressure.

Figure 8 shows the pressure solution for continuous approximation and different mesh sizes. The exact jump is ${\textstyle \Delta p=\gamma /R=4}$. We observe that the pressure solution oscillates at the interface, and it improves with finer meshes. In the case of discontinuous pressure approximation shown in Figure 9, the solution is already excellent for coarse meshes. The velocity error (measured in the norm ${\textstyle ||u||_{\infty }=\max _{i}|u_{i}|}$) for both approximations is shown in Table 1. The discontinuous pressure produces velocity solutions three orders of magnitude better than the continuous one.

(a)	(b)
Figure 8: Continuous pressure approximation: (a) profile at final time for different $h$, and (b) pressure field for $h=1/20$.

(a)	(b)
Figure 9: Discontinuous pressure approximation: (a) profile at final time for different $h$, and (b) pressure field for $h=1/20$.

Table. 1 Velocity error $||u||_{\infty }$ at final time for continuous and discontinuous pressure approximations.

${\textstyle h}$	Continuous	Discontinuous
1/20	$4.4\times {10}^{-2}$	$2.8\times {10}^{-5}$
1/40	$1.7\times {10}^{-2}$	$1.3\times {10}^{-5}$
1/80	$7.2\times {10}^{-3}$	$8.9\times {10}^{-6}$

Finally, Figure 10 shows the influence of the parameter ${\textstyle \lambda }$ on the pressure solution at the interface. For the minimum value ${\textstyle \alpha _{min}}$ that makes the system Eq. (23) invertible, the pressure profile is flat at the interface, i.e. the jump is represented exactly and incompressibility is satisfied. The larger the value of ${\textstyle \alpha }$, the more strongly the pressure continuity condition at the interface is imposed.

Figure 10: Pressure profile for discontinuous approximation and $h=1/20$. Influence of $\alpha$ value on the pressure jump at the interface: $\alpha _{min}$, $10\alpha _{min}$, $20\alpha _{min}$.

The numerical scheme developed in the previous sections will be tested in the problem of a bubble rising in a liquid column.

7 NUMERICAL EXAMPLE: RISING BUBBLE

Bubbles and bubbly flows play a significant role in a wide range of geophysical and industrial processes, such as mixing in chemical reactors, elaboration of alloys, cooling of nuclear reactors, two-phase heat exchangers, aeration processes, and atmosphere-ocean exchanges.

The shape and rise velocity of a bubble are controlled by the physical properties of the fluids and the surrounding flow field. Grace [39] developed a well-known graphical correlation for single gas bubbles rising in an infinite liquid using three dimensionless numbers: the Reynolds number (${\textstyle Re}$), the Eötvös number (${\textstyle Eo}$), and the Morton number (${\textstyle Mo}$).

$Re=\displaystyle {\frac {\rho \,L\,U}{\mu }}$ (26) $Eo=\displaystyle {\frac {\rho \,g\,L^{2}}{\gamma }}$ (27) $Mo=\displaystyle {\frac {g\,\mu ^{4}}{\rho \,\gamma ^{3}}}$ (28)

where ${\textstyle L}$ is an equivalent diameter of the bubble, and ${\textstyle U}$ the rising speed of the bubble. [39] classifies the bubble shapes into three regimes: spherical, ellipsoidal, and spherical cap. Bubbles with low ${\textstyle Re}$ or low ${\textstyle Eo}$ are spherical. For higher ${\textstyle Re}$, bubbles have an ellipsoidal shape at intermediate ${\textstyle Eo}$ and spherical cap shape at high ${\textstyle Eo}$. A more detailed regime diagram was given by Clift et al. [40], which included wobbling bubbles for ${\textstyle Re\sim 10^{3}}$ in the ellipsoidal regime, and the spherical cap regime was subdivided into spherical cap, skirted, and dimpled ellipsoidal cap for high, intermediate, and low ${\textstyle Re}$ numbers, respectively. These diagrams were further developed by Bhaga et al. [41], who also studied the flow field around the bubble, specially the wake that forms in the rear of the bubble for intermediate and high ${\textstyle Re}$.

Numerous experimental studies have been performed to understand the flow dynamics of a single rising bubble, e.g. by [42,39,43,40,41,44,45,46], but the fact that approximate theoretical solutions have only been derived in the limit of very small bubble deformations, together with the difficulties in experimentally measuring the flow variables of the bubble without any interference while it is rising and deforming, make numerical simulation an important tool to gain insight of the flow.

Previous numerical studies have mostly followed the fixed mesh approach: the pioneering works [47], [48], and [49] used the front-tracking method (together with finite differences), also [50] (finite differences) and [51] (finite volume method); level set is used by [27,26,1] (finite elements) and [52] (finite volumes); Volume-of-Fluid is used in [53] and [54]; [55] use a hybrid approach between VOF and level set (finite volumes); and interface fitting method in [56] and [45]. The Lattice-Boltzmann method is used e.g. in [57] and [58]. In these numerical works, results are qualitatively compared with experimental observations of bubble shape under different regimes. Only recently, quantitative tests for two-dimensional rising bubbles have been proposed by Hysing et al. [1]. In Section 7.1 we compare the PFEM results with these test solutions, and in Section 7.2 we propose two problems on bubble breakup and coalescence to assess the ability of a multi-fluid code to model topology changes.

7.1 Comparison with previous numerical experiments

Figure 11: Initial configuration.

The problem consists in a bubble rising in a liquid column as illustrated in Figure 11. Two tests have been proposed in [1]: test 1 considers a bubble in the ellipsoidal regime which undergoes moderate shape deformation, while in the test 2 the bubble belongs to the skirted regime and experiences much larger deformation. Both fluids are Newtonian, incompressible and isothermal, and their physical properties are listed in Table 2.

Table. 2 Physical parameters.

Test case	$\rho _{1}$	$\rho _{2}$	$\mu _{1}$	$\mu _{2}$	g	$\gamma$	$Re$	$Eo$	$\rho _{1}/\rho _{2}$	$\mu _{1}/\mu _{2}$
1	1000	100	10	1	0.98	24.5	35	10	10	10
2	1000	1	10	0.1	0.98	1.96	35	125	1000	100

The reference solutions presented in [1] have been run with three different numerical approaches: the TP2D [59], the FreeLIFE [60], and the MooNMD [61]. They all use the finite element method, but the two first approaches describe the interface with the level set, while the latter tracks it in an arbitrary Lagrangian-Eulerian way. More specific details about the methods can be found in the original publication. The following bubble quantities are used to compare the results:

• shape at the final time ${\textstyle t=3}$ s,
• circularity ${\textstyle \displaystyle c={\frac {2{\sqrt {\pi {\mbox{Area}}}}}{\mbox{Perimeter}}}}$,
• center of mass ${\textstyle \displaystyle {\boldsymbol {X}}_{c}=\int _{\Omega _{2}}{\boldsymbol {x}}\,dx/\int _{\Omega _{2}}1\,dx}$, and
• rise velocity ${\textstyle \displaystyle {\boldsymbol {U}}_{c}=\int _{\Omega _{2}}u\,dx/\int _{\Omega _{2}}1\,dx}$.

The computations have been performed on unstructured meshes with element size ${\textstyle h=1/40}$ in the bulk of the fluids and wall regions, and the mesh at the interface has been refined to ${\textstyle h=1/80}$, ${\textstyle 1/160}$, ${\textstyle 1/320}$ and ${\textstyle 1/640}$ in order to analyze the convergence in ${\textstyle h}$ of the solution. With a refinement based on the distance to the interface (see Figure 12), we can use an arbitrarily fine mesh without increasing the total number of nodes to impractical values as would be the case with an uniform mesh.

Figure 12: Mesh is refined close to the interface with the help of a distance function.

(a) ${\textstyle t=0}$ s (b) ${\textstyle t=0.5}$ s (c) ${\textstyle t=1.0}$ s (d) ${\textstyle t=1.5}$ s (e) ${\textstyle t=2.0}$ s (f) ${\textstyle t=2.5}$ s (g) ${\textstyle t=3.0}$ s

Figure 13: Test 1. Bubble evolution for mesh size $h=1/320$.

For test 1, Figure 13 shows the evolution of the rising bubble. At final time ${\textstyle t=3}$ s we compare the PFEM bubble shapes for the meshes ${\textstyle h=1/40,1/80,1/160}$ and ${\textstyle 1/320}$, and observe that they converge to the shape of the finest mesh (Figure 14), which is in good agreement with the TP2D solution reported in [1] (Figure 15a). The plots of the bubble circularity (Figure 15b) and rise velocity (Figure 15d) show that our bubble is slightly oscillating, but the evolution of the center of mass (Figure 15c) is again in good agreement. The oscillating behavior of the PFEM results may be explained by the fact that, on the one side, PFEM does not introduce diffusivity at the interface, and on the other side, the geometrical method we use to calculate the curvature (the osculating circle) may not be accurate enough. Regarding the volume conservation, without any correction technique the bubble volume variation between the initial and final times, ${\textstyle \Delta V={\frac {V_{f}-V_{0}}{V_{0}}}}$, is of order ${\textstyle {\cal {O}}(10^{-4})}$ (Table 3).

Figure 14: Test 1. PFEM bubble shape for different mesh sizes at $t=3$ s.

(a) Shape at ${\textstyle t=3}$ s	(b) Circularity
(c) Center of mass	(d) Rise velocity
Figure 15: Test 1. Comparison of benchmark quantities: PFEM ($h=1/320$) vs. TP2D results.

Table. 3 Volume variation

Test	$h$	Initial Volume	Final Volume	Variation
1	1/320	0.19635	0.1965	$+7\times {10}^{-4}$
2	1/640	0.19635	0.19633	$-1\times {10}^{-4}$

Same type of results are shown for test 2 in Figures 16, 17, and 18. Although the bubbles in both test cases rise with similar velocity, the decrease in surface tension causes bubble 2 to undergo a much larger deformation and to develop thin filaments. In the TP2D and FreeLIFE solutions these filaments break up, in contrast to the moving mesh solutions of PFEM and MooNMD (Figure 18a). In the physical reality, breakup occurs due to capillary waves present on the interface, which trigger the three-dimensional Plateau-Rayleigh instability when the filament radius is small enough. Thus, capillary waves can cause the skirt filament to fragment during flow, though this response requires very large elongations, typically greater than 20 times the initial bubble radius [62]. Figure 19 shows the mesh of the PFEM solution in the skirted region. The filament is not thin enough to break up. The volume variation is excellent again, of order ${\textstyle {\cal {O}}(10^{-4})}$.

(16) Test 2. Bubble evolution for mesh size $h=1/640$.

(a) ${\textstyle t=0}$ s (b) ${\textstyle t=0.5}$ s (c) ${\textstyle t=1.0}$ s (d) ${\textstyle t=1.5}$ s (e) ${\textstyle t=2.0}$ s (f) ${\textstyle t=2.5}$ s (g) ${\textstyle t=3.0}$ s

Figure 17: Test 2. PFEM bubble shape for different mesh sizes at $t=3$ s.

(a) Shape at ${\textstyle t=3}$ s	(b) Circularity
(c) Center of mass	(d) Rise velocity
Figure 18: Test 2. Benchmark quantities comparison of PFEM (black line) and TP2D (red), FreeLIFE (green) and MooNMD (blue) results.

Figure 19: Test 2. Detail of bubble skirt at $t=3$ s.

The fact of using an explicit approach when discretizing in time the surface tension force introduces the stability condition for the time step expressed in Eq. (24). If larger time steps are taken, instabilities will develop at the interface. This time step constraint is very restrictive for fine meshes (${\textstyle \Delta t_{st}\propto h^{3/2}}$). In our case, the refined mesh close to the interface imposes rather small global time steps (for test 1 with mesh ${\textstyle h=1/320}$, ${\textstyle \Delta t_{st}<3.3\times 10^{-4}}$; and for test 2 with ${\textstyle h=1/640}$, ${\textstyle \Delta t_{st}<3.9\times 10^{-4}}$) that undoubtedly affect the computational efficiency.

7.2 Bubble breakup and coalescence

Topology changes in multi-fluid flows can be divided into two classes [63]:

• Films that fragment. If a bubble approaches another bubble or a flat surface, the fluid in between must be squeezed out before the bubbles are sufficiently close so that the film becomes unstable to attractive forces and fragment.
• Threads that break. A long and thin cylinder of fluid will generally break by the Plateau-Rayleigh instability in the region where the cylinder becomes sufficiently thin so that surface tension pinches it into two.

In order to test the capabilities of PFEM to handle interfaces with changing topology, and motivated by the disagreement in the solution of test 2, we have simulated two examples on a film that fragments, namely the breakup and coalescence of bubbles.

We consider the same fluid properties and configuration of test 1. In the case of the breakup, we add a flat interface at ${\textstyle y=1}$ so that the upper region belongs to the same fluid than the bubble (see Figure 20a). The bubble rises, approaching the flat interface. The film of heavy fluid that separates the two regions of light fluid becomes thinner and thinner until it fragments and the regions fuse (Figure 20). Whereas in the physical reality the fragmentation of the film is caused by attractive forces at the microscopic scale (forces which are usually not included in the continuum description), in our simulations fragmentation is caused by a connectivity change at the interface, as illustrated in Figure 21.

(a) ${\textstyle t=0}$ s (b) ${\textstyle t=1.5}$ s (c) ${\textstyle t=2.5}$ s (d) ${\textstyle t=3.5}$ s (e) ${\textstyle t=4.5}$ s (f) ${\textstyle t=5.5}$ s (g) ${\textstyle t=6.0}$ s (h) ${\textstyle t=6.5}$ s

Figure 20: Bubble breakup.

One of the main difficulties we face in our Lagrangian approach is the connectivity changes introduced by the remeshing process. In general, these reconnections may alter the equilibrium at the interface, slow down convergence and affect mass conservation. Thus, in interfacial flows it is essential to avoid them. We are using an unconstrained Delaunay triangulator which does not allow to fix connectivities. Therefore, to ensure that a specific connectivity remains, we refine long interfacial edges and remove nodes too close to the interface. Unfortunately, this strategy would preclude the possibility of breakup, as the interface could elongate endlessly. In the way PFEM defines interfaces, it is possible to have fluid regions spanned by just one element layer. The breakup criterium we have implemented in PFEM is to permit connectivity changes in elements where all nodes lie at the interface. In this way, a thin fluid thread can stop elongating and fragment. Breakup is then dependent on the mesh resolution, that is, it happens when the thickness of the film is similar to the mesh resolution of the interface. This is not a drawback specific of PFEM, breakup is mesh dependent in front-capturing methods as well. For example, in the level set method, two interfaces are described as two different zero contours of the same level set function, and these interfaces will automatically merge once they get close enough, relative to the spatial resolution of the mesh where the level set function is defined. The resolution determines the smallest distance between two zero level sets of the level set function for which they can still be distinguished as separate zero contours. Interfaces can in fact merge faster due to the diffusivity of the schemes used for advection and reinitialization of the level set function [27].

(a) ${\textstyle t^{n}}$	(b) ${\textstyle t^{n+1}}$
Figure 21: Connectivity change that produces breakup at fluid films spanned by just one mesh element.

The pressure field at and after breakup is shown in Figures 22 and 23. The different pressure values inside and outside the bubble equilibrate after breakup, what occurs at ${\textstyle t=5.97}$ s.

(a) ${\textstyle t=5.965}$ s	(b) ${\textstyle t=5.97}$ s
(c) ${\textstyle t=5.975}$ s	(d) ${\textstyle t=5.98}$ s
Figure 22: Pressure field at breakup (variable scale ranges in legend).

(a) ${\textstyle t=6.0}$ s	(b) ${\textstyle t=6.05}$ s
(c) ${\textstyle t=6.125}$ s	(d) ${\textstyle t=6.18}$ s
Figure 23: Pressure field after breakup (variable scale ranges in legend).

For the simulation of bubble coalescence, we consider the same rectangular domain ${\textstyle (0,1)\times (0,2)}$ as before, with two circular bubbles inside. The center of the first bubble is ${\textstyle (0.5,1.0)}$ and its radius is equal to 0.25, the center of the second bubble is ${\textstyle (0.5,0.5)}$ and the radius 0.2. Since the small bubble is located close to the large one, this lower bubble turns out to be in the wake of the upper bubble and rises faster than that. Figure 24 shows the coalescence process. The mechanism is again the rupture of the thin film between the bubbles. This happens not during the impact of the bubbles (around ${\textstyle t=2.5}$ s) but during the separation after impact, as corresponds to the physical reality [64].

(a) ${\textstyle t=0}$ s (b) ${\textstyle t=0.5}$ s (c) ${\textstyle t=1.0}$ s (d) ${\textstyle t=1.5}$ s (e) ${\textstyle t=2.0}$ s (f) ${\textstyle t=2.5}$ s (g) ${\textstyle t=3.0}$ s (h) ${\textstyle t=3.5}$ s (i) ${\textstyle t=4.0}$ s (j) ${\textstyle t=4.5}$ s (k) ${\textstyle t=5.0}$ s (l) ${\textstyle t=10.0}$ s

Figure 24: Bubble coalescence.

8 SUMMARY AND CONCLUSIONS

We have developed a numerical scheme for the simulation of multi-fluid flows with the Particle Finite Element Method able to deal with large jumps in the physical properties (included surface tension), and able to accurately represent the ${\textstyle {\cal {C}}^{-1}}$ and ${\textstyle {\cal {C}}^{0}}$ discontinuities in the flow variables. Interfaces are tracked by the moving mesh, pressure degrees of freedom have been duplicated at the interface nodes to represent the ${\textstyle {\cal {C}}^{-1}}$ discontinuity of this variable due to surface tension and variable viscosity (Eq. 9), velocity and pressure variables have been decoupled through a pressure segregation method, and the mesh has been refined in the vicinity of the interface to improve the accuracy and efficiency of the computations.

We have first tested the scheme in a simple static bubble problem and compared the results for continuous and discontinuous pressure approximations. We have then applied the method to the more complicated rising bubble problem presented in Hysing et al. [1]. We can conclude that PFEM solutions for the single rising bubble are in good agreement with those reported in [1]. For test 1, our bubble is slightly oscillating in contrast to the reference solution. A reason for this may be that PFEM does not introduce diffusion at the interface. In any case, more comparisons with other methods are needed. For test 2, although PFEM can handle interface breakup without problems (as shown in Section 7.2), the skirt filaments remain intact. Breakup happens only when the fluid region is spanned by just one element layer. This allows to model thin films of thickness ${\textstyle h}$, being ${\textstyle h}$ the mesh size at the interface. In both tests, we have achieved an excellent mass conservation without any correction.

Finally, we propose two bubble breakup and coalescence problems as benchmarks for testing the capabilities of multi-fluid flow codes to handle topology changes of the interface.

BIBLIOGRAPHY

[1] Hysing, S. and Turek, S. and D. Kuzmin and N. Parolini and E. Burman and S. Ganesan and L. Tobiska. (2009) "Quantitative benchmark computations of two-dimensional bubble dynamics", Volume 60. International Journal for Numerical Methods in Fluids 1259-1288

[2] Idelsohn, S.R. and Mier-Torrecilla, M. and Oñate, E. (2009) "Multi-fluid flows with the Particle Finite Element Method", Volume 198. Computer Methods in Applied Mechanics and Engineering 2750-2767

[3] G.K. Batchelor. (1967) "An introduction to fluid dynamics". Cambridge University Press

[4] Idelsohn, S.R. and Mier-Torrecilla, M. and Nigro, N. and Oñate, E. (2009) "On the analysis of heterogeneous fluids with jumps in the viscosity using a discontinuous pressure field", Volume In press. Computational Mechanics

[5] Unverdi, S.O. and Tryggvason, G. (1992) "Computations of multi-fluid flows", Volume 60. Physica D: Nonlinear Phenomena 70-83

[6] W. Shyy. (1996) "Computational Fluid Dynamics with moving boundaries". Taylor&Francis

[7] Scardovelli, R. and Zaleski, S. (1999) "Direct Numerical Simulation of free-surface and interfacial flow", Volume 31. Annual Reviews of Fluid Mechanics 567-603

[8] Caboussat, A. (2005) "Numerical simulation of two-phase free surface flows", Volume 12. Archives of Computational Methods in Engineering 165-224

[9] Hughes, T.J.R. and Liu, W.K. and Zimmermann, T.K. (1981) "Lagrangian-Eulerian finite element element formulation for incompressible viscous flows", Volume 29. Computer Methods in Applied Mechanics and Engineering 239-349

[10] B. Ramaswamy and M. Kawahara. (1986) "Arbitrary Lagrangian-Eulerian finite element method for the analysis of free surface fluid flows", Volume 1. Computational Mechanics 103-108

[11] T.E. Tezduyar and M. Behr and J. Liou. (1992) "A new strategy for finite element computations involving moving boundaries and interfaces. The deforming-spatial-domain/space-time procedure: I. The concept and preliminary numerical tests", Volume 94. Computer Methods in Applied Mechanics and Engineering 339-351

[12] T.E. Tezduyar and M. Behr and S. Mittal and J. Liou. (1992) "A new strategy for finite element computations involving moving boundaries and interfaces. The deforming-spatial-domain/space-time procedure: II. Computation of free-surface flows, two-liquid flows and flows with drifting cylinders", Volume 94. Computer Methods in Applied Mechanics and Engineering 353-371

[13] Hirt, C.W. and Cook, J.L. and Butler, T.D. (1970) "A Lagrangian method for calculating the dynamics of an incompressible fluid with free surface", Volume 5. Journal of Computational Physics 103-124

[14] B. Ramaswamy and M. Kawahara. (1987) "Lagrangian finite element analysis applied to viscous free surface fluid flow", Volume 7. International Journal for Numerical Methods in Fluids 953-984

[15] Idelsohn, S.R. and Oñate, E. and Del Pin, F. (2003) "A Lagrangian meshless finite element method applied to fluid-structure interaction problems", Volume 81. Computers and Structures 8–11 655–671

[16] Idelsohn, S.R. and Oñate, E. and Del Pin, F. (2004) "The Particle Finite Element Method: A powerful tool to solve incompressible flows with free-surfaces and breaking waves", Volume 61. International Journal for Numerical Methods in Engineering 7 964–989

[17] Oñate, E. and Idelsohn, S.R. and Del Pin, F. and Aubry, R. (2004) "The Particle Finite Element Method: An Overview", Volume 1. International Journal of Computational Methods 2 267–307

[18] Del Pin, F. and Idelsohn, S.R. and Oñate, E. and Aubry, R. (2007) "The ALE/Lagrangian Particle Finite Element Method: A new approach to computation of free-surfaces flows and fluid-object interactions", Volume 36. Computers and Fluids 1 27–38

[19] Harlow, F.H. (1964) "The Particle-in-Cell computing method for fluid dynamics", Volume 3. Methods in Computational Physics 313-343

[20] Unverdi, S.O. and Tryggvason, G. (1992) "A front-tracking method for viscous, incompressible, multi-fluid flows", Volume 100. Journal of Computational Physics 25-37

[21] Hirt, C. and Nichols, B. (1981) "Volume of Fluid (VOF) method for the dynamics of free boundaries", Volume 39. Journal of Computational Physics 201-225

[22] Osher, S. and Sethian, J. (1988) "Fronts propagating with curvature dependant speed: algorithms based on Hamilton-Jacobi formulations", Volume 79. Journal of Computational Physics 12-49

[23] Idelsohn, S.R. and Calvo, N. and Oñate, E. (2003) "Polyhedrization of an arbitrary 3D point set", Volume 192. Computer Methods in Applied Mechanics and Engineering 2649-2667

[24] Oñate, E. and Idelsohn, S.R. and Celigueta, M.A. and Rossi, R. (2008) "Advances in the Particle Finite Element Method for the analysis of fluid-multibody interaction and bed erosion in free surface flows", Volume 197. Computer Methods in Applied Mechanics and Engineering 1777-1800

[25] van Kan, J. (1986) "A second-order accurate pressure correction scheme for viscous incompressible flow", Volume 7. SIAM Journal on Scientific and Statistical Computing 870-891

[26] A. Smolianski. (2001) "Numerical Modeling of Two Fluid Interfacial Flows". University of Jyväskylä

[27] A.-K. Tornberg. (2000) "Interface Tracking Methods with Application to Multiphase Flows". Royal Institute of Technology, Stockholm

[28] Minev, P.D. and Chen, T. and Nandakumar, K. (2003) "A finite element technique for multifluid incompressible flow using Eulerian grids", Volume 187. Journal of Computational Physics 255-273

[29] Coppola-Owen, A.H. and Codina, R. (2005) "Improving Eulerian two-phase flow finite element approximation with discontinuous gradient pressure shape functions", Volume 49. International Journal for Numerical Methods in Fluids 1287-1304

[30] J. Chessa and T. Belytschko. (2003) "An extended finite element method for two-phase fluids", Volume 70. ASME Journal of Applied Mechanics 10-17

[31] S. Gross and A. Reusken. (2007) "An extended pressure finite element space for two-phase incompressible flows with surface tension", Volume 224. Journal of Computational Physics 40-58

[32] Juntunen, M. and Stenberg, R. (2007) "Nitsche's method for general boundary conditions". Research Reports A530, University of Technology, Institute of Mathematics

[33] Mier-Torrecilla, M. (2010) "Numerical simulation of multi-fluid flows with the Particle Finite Element Method". Technical University of Catalonia

[34] J.U. Brackbill and D.B. Kothe and C. Zemach. (1992) "A continuum method for modeling surface tension", Volume 100. Journal of Computational Physics 335-354

[35] Bänsch, E. (2001) "Finite Element discretization of the Navier-Stokes equations with a free capillary surface", Volume 88. Numerische Mathematik 203-235

[36] Hysing, S. (2006) "A new implicit surface tension implementation for interfacial flows", Volume 51. International Journal for Numerical Methods in Fluids 659-672

[37] Popinet, S. and Zaleski, S. (1999) "A front tracking algorithm for accurate representation of surface tension", Volume 30. International Journal for Numerical Methods in Fluids 775-793

[38] Sussman, M. and Smereka, P. and Osher, S. (1994) "A level set approach for computing solutions to incompressible two-phase flows", Volume 114. Journal of Computational Physics 146-159

[39] Grace, J. R. (1973) "Shapes and Velocities of Bubbles Rising in Infinite Liquids", Volume 51. Transactions of the Institution of Chemical Engineers 116-120

[40] R. Clift and J.R. Grace and M.E. Weber. (1978) "Bubbles, drops and particles". Academic Press

[41] Bhaga, D. and Weber, M.E. (1981) "Bubbles in viscous liquids: shapes, wakes and velocities", Volume 105. Journal of Fluid Mechanics 61-85

[42] Haberman, W. L. and Morton, R. K. (1954) "An experimental study of bubbles moving in liquids", Volume 387. Proceedings of the American Society of Civil Engineers 1-25

[43] Hnat, J. G. and Buckmaster, J. D. (1976) "Spherical cap bubbles and skirt formation", Volume 19. Physics of Fluids 182-194

[44] Maxworthy, T. and Gnann, C. and Kuerten, M. and Durst, F. (1996) "Experiments on the rise of air bubbles in clean viscous liquids", Volume 321. Journal of Fluid Mechanics 421-441

[45] F. Raymond and J.M. Rosant. (2000) "A numerical and experimental study of the terminal velocity and shape of bubbles in viscous liquids", Volume 55. Chemical Engineering Science 943-955

[46] M. Wu and M. Gharib. (2002) "Experimental Studies on the Shape and Path of Small Air Bubbles Rising in Clean Water", Volume 14. Physics of Fluids 49-52

[47] Esmaeeli, A. and Tryggvason, G. (1998) "Direct numerical simulations of bubbly flows. Part 1. Low Reynolds number arrays", Volume 377. Journal of Fluid Mechanics 313-345

[48] Esmaeeli, A. and Tryggvason, G. (1999) "Direct numerical simulations of bubbly flows. Part 2. Moderate Reynolds number arrays", Volume 385. Journal of Fluid Mechanics 325-358

[49] Bunner, B. and Tryggvason, G. (2002) "Dynamics of homogeneous bubbly flows. Part 1. Rise velocity and microstructure of the bubbles", Volume 466. Journal of Fluid Mechanics 17-52

[50] F.S. de Sousa and N. Mangiavacchi and L.G. Nonato and A. Castelo and M.F. Tomé and V.G. Ferreira and J.A. Cuminato and S. McKee. (2004) "A front-tracking/front-capturing method for the simulation of 3D multi-fluid flows with free surfaces", Volume 198. Journal of Computational Physics 469-499

[51] Hua, J. and Lou, J. (2007) "Numerical simulation of bubble rising in viscous liquid", Volume 222. Journal of Computational Physics 769-795

[52] Z. Yu and L.-S. Fan. (2008) "Direct Simulation of the Buoyant Rise of Bubbles in Infinite Liquid Using Level Set Method", Volume 86. Canadian Journal of Chemical Engineering 267-275

[53] L. Chen and S.V. Garimella and J.A. Reizes and E. Leonardi. (1999) "The development of a bubble rising in a viscous liquid", Volume 387. Journal of Fluid Mechanics 61-96

[54] Van Sint Annaland, M. and N. G. Deen and J. A. M. Kuipers. (2005) "Numerical Simulation of Gas Bubbles Behavior Using a Three-Dimensional Volume of Fluid Method", Volume 60. Chemical Engineering Science 2999-3011

[55] Bonometti, T. and Magnaudet, J. (2007) "An interface-capturing method for incompressible two-phase flows. Validation and application to bubble dynamics", Volume 33. International Journal of Multiphase Flow 109-133

[56] Ryskin, G. and Leal, L. G. (1984) "Numerical solution of free-boundary problems in fluid mechanics. Part 2. Buoyancy-driven motion of a gas bubble through a quiescent liquid", Volume 148. Journal of Fluid Mechanics 19-35

[57] X. Frank and D. Funfschilling and N. Midoux and H. Z. Li. (2006) "Bubbles in a viscous liquid: Lattice Boltzmann simulation and experimental validation", Volume 546. Journal of Fluid Mechanics 113-122

[58] I. O. Kurtoglu and C. L. Lin. (2006) "Lattice Boltzmann Study of Bubble Dynamics", Volume 50. Numerical Heat Transfer B 333-351

[59] S. Turek. (1998) "Efficient Solvers for Incompressible Flow Problems". Springer

[60] Parolini, N. and Burman, E. (2005) "A finite element level set method for viscous free-surface flows". 7th Conference on Applied and Industrial Mathematics in Italy 416-427

[61] Ganesan, S. and Matthies, G. and Tobiska, L. (2007) "On spurious velocities in incompressible flow problem with interfaces", Volume 196. Computer Methods in Applied Mechanics and Engineering 1193-1202

[62] Stone, H. (1994) "Dynamics of drop deformation and breakup in viscous fluids", Volume 26. Annual Reviews of Fluid Mechanics 65-102

[63] Tryggvason, G. and Bunner, B. and Esmaeeli, A. and Juric, D. and Al-Rawahi, N. and Tauber, W. and Han, J. and Nas, S. and Jan, Y.-J. (2001) "A front-tracking method for the computations of multiphase flow", Volume 169. Journal of Computational Physics 708-759

[64] N. Bremond and A. Thiam and J. Bibette. (2008) "Decompressing Emulsion Droplets Favors Coalescence", Volume 100. Physical Review Letters