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

​

'''M.&nbsp;Mier-Torrecilla, S.R.&nbsp;Idelsohn<span id="fnc-1"></span>[[#fn-1|<sup>1</sup>]]and E.&nbsp;Oñate'''

-->

​

==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 <span id='citeF-1'></span>[[#cite-1|[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:

​

<ol>

​

<li>Accurate definition of the interface position. 

​

The interface separating the fluids needs to be tracked accurately without introducing excessive numerical smoothing. </li>

<li>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. </li>

<li>Modeling the discontinuities of the flow variables across the interface.

​

Velocity and pressure may be discontinuous across the interface under certain conditions. </li>

<li>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. </li>

​

</ol>

​

This paper extends the work of the authors in the study of multi-fluid flows with the Particle Finite Element Method <span id='citeF-2'></span>[[#cite-2|[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 GOVERNING EQUATIONS|2]] we present the governing equations together with the interface conditions and the possible discontinuities of the flow variables. In Section [[#3 INTERFACE DESCRIPTION|3]] we briefly review the interface descriptions most used in the literature, to focus later in the Particle Finite Element Method (Section [[#4 PARTICLE FINITE ELEMENT METHOD|4]]) and the numerical scheme we have developed (Sections [[#5 PRESSURE SEGREGATION FOR THE MULTI-FLUID EQUATIONS|5]] and [[#6 SURFACE TENSION|6]]). Finally, in Section [[#7 NUMERICAL EXAMPLE: RISING BUBBLE|7]] we apply the PFEM to the solution of the two-dimensional bubble rise, breakup and coalescence problems.

​

==2 GOVERNING EQUATIONS==

​

<div id='img-1'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-NSdomain.png|180px|Two-fluid flow configuration.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 1:''' Two-fluid flow configuration.

|}

​

Let <math display="inline">\Omega \subset \boldsymbol{R}^d</math>, <math display="inline">d\in \{ 2,3\} </math>, be a bounded domain containing two different fluids (see Figure [[#img-1|1]]). We denote time by <math display="inline">t</math>, the Cartesian spatial coordinates by <math display="inline">\boldsymbol{x}=\{ x_i\} _{i=1}^d</math>, and the vectorial operator of spatial derivatives by <math display="inline">\boldsymbol{\nabla }=\{ \partial _{x_i}\} _{i=1}^d</math>. The evolution of the velocity <math display="inline">u=u(\boldsymbol{x},t)</math> and the pressure <math display="inline">p=p(\boldsymbol{x},t)</math> is governed by the Navier-Stokes equations:

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\rho \displaystyle \frac{{d} u}{{d} t}=\boldsymbol{\nabla }\cdot{\boldsymbol{\sigma }}+\rho \boldsymbol{g}\quad \mbox{ in} \;\Omega \times (0,T) </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (1.a)

|-

| style="text-align: center;" | <math> \displaystyle \frac{{d} \rho }{{d} t}+\rho \boldsymbol{\nabla }\cdot{u}=0  \quad \mbox{ in} \;\Omega \times (0,T) </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (1.b)

|}

|}

​

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

​

<span id="eq-2"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\boldsymbol{\sigma }=-p\boldsymbol{I}+2\mu (\boldsymbol{D}-\frac{1}{3}\varepsilon _V\boldsymbol{I})  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (2)

|}

​

with <math display="inline">\boldsymbol{I}</math> the identity tensor, <math display="inline">\mu </math> the dynamic viscosity, <math display="inline">\boldsymbol{D}=\frac{1}{2}(\boldsymbol{\nabla }u+\boldsymbol{\nabla }^Tu)</math> the strain rate tensor, and <math display="inline">\varepsilon _V=\boldsymbol{\nabla }\cdot{u}</math> the volumetric strain rate.

​

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

​

<span id="eq-3"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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.  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (3)

|}

​

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

​

===Boundary and interface conditions===

​

In order for the Navier-Stokes problem ([[#2 GOVERNING EQUATIONS|2]]) to be well-posed, suitable boundary conditions need to be specified. On the external boundary <math display="inline">\partial \Omega =\Gamma _D\cup \Gamma _N</math>, such that <math display="inline">\Gamma _D\cap \Gamma _N=\emptyset </math>, we consider the following:

​

<span id="eq-4"></span>

<span id="eq-5"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>u=\bar u  \quad \mbox{ on} \;\Gamma _D </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (4)

|-

| style="text-align: center;" | <math> \boldsymbol{\sigma }\cdot \boldsymbol{n}=\bar \boldsymbol{\sigma }_n   \quad \mbox{ on} \;\Gamma _N   </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (5)

|}

|}

​

<math display="inline">\bar u</math> is the prescribed velocity, <math display="inline">\boldsymbol{n}</math> the outer unit normal to <math display="inline">\Gamma _N</math>, and <math display="inline">\bar \boldsymbol{\sigma }_n</math> the prescribed traction vector. A Neumann boundary <math display="inline">\Gamma _N</math> with <math display="inline">\bar \boldsymbol{\sigma }_n=\boldsymbol{0}</math> is called ''free surface''.

​

On the internal interfaces <math display="inline">\Gamma _{int}</math>, the coupling conditions are <span id='citeF-3'></span>[[#cite-3|[3]]]:

​

<span id="eq-6"></span>

<span id="eq-7"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\mbox{⟦}u\mbox{⟧}=\boldsymbol{0}  \qquad \mbox{ on} \;\Gamma _{int} </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (6)

|-

| style="text-align: center;" | <math> \mbox{⟦}\boldsymbol{\sigma }\mbox{⟧}\cdot \boldsymbol{n}=\gamma \kappa \boldsymbol{n}  \qquad \mbox{ on} \;\Gamma _{int}  </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (7)

|}

|}

​

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

​

Equation ([[#eq-6|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 (<math display="inline">\mu ^+,\mu ^->0</math>), similar to a no-slip condition.

​

Equation ([[#eq-7|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 <math display="inline">\Gamma _{int}</math>. The surface tension coefficient <math display="inline">\gamma </math> 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:

​

* ''<math>{\cal C}^0</math> discontinuity'', when the flow variable has a kink (i.e.&nbsp;the gradient has a jump), and

* ''<math>{\cal C}^{-1}</math> 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 <math display="inline">{\cal C}^0</math> discontinuity in the pressure field (Fig.&nbsp;[[#img-2|2]]a).

​

Differences in viscosity lead to discontinuous components of the strain rate tensor <math display="inline">\boldsymbol{D}</math>, and therefore to a <math display="inline">{\cal C}^0</math> discontinuity of the velocity field at the interface (Fig.&nbsp;[[#img-2|2]]b):

​

<span id="eq-8"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (8)

|}

​

with <math display="inline">\partial _s = \boldsymbol{t}\cdot \boldsymbol{\nabla }</math> the tangential derivative.

​

Both differences in viscosity and the presence of surface tension cause a <math display="inline">{\cal C}^{-1}</math> discontinuity in the pressure field (Figs.&nbsp;[[#img-2|2]]b and [[#img-2|2]]c), as shown in <span id='citeF-4'></span>[[#cite-4|[4]]]:

​

<span id="eq-9"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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   </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (9)

|}

​

Notice that even in the case of <math display="inline">\gamma=0</math>, pressure is discontinuous when <math display="inline">\mu ^+\neq \mu ^-</math>.

​

<div id='img-2a'></div>

<div id='img-2b'></div>

<div id='img-2c'></div>

<div id='img-2'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-discontinuities_density.png|372px|(a)]]

|[[Image:draft_Samper_466629376-discontinuities_viscosity.png|372px|(b)]]

|- style="text-align: center; font-size: 75%;"

| (a) 

| (b) 

|-

| colspan="2"|[[Image:draft_Samper_466629376-discontinuities_surfacetension.png|372px|(c)]]

|- style="text-align: center; font-size: 75%;"

|  colspan="2" | (c) 

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''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.&nbsp;<span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-7'></span><span id='citeF-8'></span>[[#cite-5|[5,6,7,8]]]. <div id='img-3a'></div>

<div id='img-3b'></div>

<div id='img-3'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-interfacemethods_moving.png|180px|(a)]]

|[[Image:draft_Samper_466629376-interfacemethods_fixed.png|180px|(b)]]

|- style="text-align: center; font-size: 75%;"

| (a) 

| (b) 

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''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 [[#img-3|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 <span id='citeF-9'></span><span id='citeF-10'></span>[[#cite-9|[9,10]]], the deformable-spatial-domain/stabilized space-time deformation (DSD/SST) method <span id='citeF-11'></span><span id='citeF-12'></span>[[#cite-11|[11,12]]], or the fully Lagrangian formulation such as in <span id='citeF-13'></span><span id='citeF-14'></span>[[#cite-13|[13,14]]] and the Particle Finite Element Method <span id='citeF-15'></span><span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-18'></span>[[#cite-15|[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 <span id='citeF-19'></span><span id='citeF-20'></span>[[#cite-19|[19,20]]], and ''front-capturing methods'', which introduce a new variable <math display="inline">\psi </math> 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 <span id='citeF-21'></span>[[#cite-21|[21]]], and the Level Set method by Osher et al.&nbsp;<span id='citeF-22'></span>[[#cite-22|[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.&nbsp;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 <span id='citeF-23'></span>[[#cite-23|[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.&nbsp;<span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-24'></span>[[#cite-16|[16,17,24]]].

​

Since in PFEM the interface is described by mesh nodes and element edges (Fig. [[#img-3|3]]a), 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.

​

<div id='img-4'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-pressureprofiles.png|420px|Pressure profiles when using continuous and discontinuous representations.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''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 Interface discontinuities|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, <math display="inline">{\cal C}^0</math> discontinuities need no special attention when the interface is aligned with the mesh, as the kinks in the solution are automatically represented.  Only <math display="inline">{\cal C}^{-1}</math> discontinuities need some attention in PFEM. In particular, the pressure field has been made double-valued at the interface, i.e.&nbsp;pressure degrees of freedom have been duplicated (<math display="inline">p^+</math>, <math display="inline">p^-</math>) in the interface nodes <span id='citeF-4'></span>[[#cite-4|[4]]]. The pressure discontinuity is thus optimally approximated. Figure [[#img-4|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)

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\displaystyle \frac{{d} \boldsymbol{x}}{{d} t} = u(\boldsymbol{x},t) </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (10)

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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 </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (11.a)

|-

| style="text-align: center;" | <math>   \boldsymbol{\nabla }\cdot{u}=0 \quad \mbox{ in}~\Omega </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (11.b)

|-

| style="text-align: center;" | <math>   u=\bar u\quad \mbox{ on}~\Gamma _D </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (11.c)

|-

| style="text-align: center;" | <math>   \boldsymbol{\sigma }\cdot \boldsymbol{n}=\bar \boldsymbol{\sigma }_n \quad \mbox{ on} \;\Gamma _N </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (11.d)

|-

| style="text-align: center;" | <math>   \mbox{⟦}u\mbox{⟧}= 0 \quad \mbox{ and} \quad \mbox{⟦}\boldsymbol{\sigma }\mbox{⟧}\cdot \boldsymbol{n}=\gamma \kappa \boldsymbol{n} \quad \mbox{ on}~\Gamma _{int} </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (11.e)

|-

| style="text-align: center;" | <math>   u(\boldsymbol{x},0) = u^0(\boldsymbol{x}) \quad \mbox{ in}~\Omega  </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (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 <span id='citeF-25'></span>[[#cite-25|[25]]], and after implicit backward Euler time discretization, we propose the following splitting of equations ([[#5 PRESSURE SEGREGATION FOR THE MULTI-FLUID EQUATIONS|5]]) <span id='citeF-26'></span>[[#cite-26|[26]]]:

​

<ol>

​

<li>Find an intermediate velocity <math display="inline">\tilde{\boldsymbol{u}}</math> solution of

​

<span id="eq-12.a"></span>

<span id="eq-12.b"></span>

<span id="eq-12.c"></span>

<span id="eq-12.d"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>

​

\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}</math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (12.a)

|-

| style="text-align: center;" | <math> \tilde{\boldsymbol{u}}=\bar u\quad \mbox{ on}~\Gamma _D </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (12.b)

|-

| style="text-align: center;" | <math>  (-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 </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (12.c)

|-

| style="text-align: center;" | <math> \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>

​

</math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (12.d)

|}

|}

<li>Determine <math display="inline">p^{n+1}</math> as solution of

​

<span id="eq-13.a"></span>

<span id="eq-13.b"></span>

<span id="eq-13.c"></span>

<span id="eq-13.d"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>

​

\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 , </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (13.a)

|-

| style="text-align: center;" | <math> \mbox{ and} \quad \boldsymbol{n}\cdot \boldsymbol{\nabla }(p^{n+1}-p^n)=0 \quad{on}~\Gamma _D  </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (13.b)

|-

| style="text-align: center;" | <math> (p^{n+1}-p^n)\boldsymbol{n}=\boldsymbol{0}\quad \mbox{ on}~\Gamma _N </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (13.c)

|-

| style="text-align: center;" | <math>  \mbox{⟦}p^{n+1}-p^n\mbox{⟧}=0 \quad \mbox{ on}~\Gamma _{int}  </li>

​

</math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (13.d)

|}

|}

​

<li>Finally, the <math display="inline">{\boldsymbol{u}}^{n+1}</math> velocity is obtained by

​

<span id="eq-14"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>

​

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

​

</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (14)

|}

​

</ol>

​

The interfacial stress jump condition has been consistently split between steps 1 and 2. Although the continuity of the pressure difference <math display="inline">p^{n+1}-p^n</math> across the interface expressed in Eq.&nbsp;([[#eq-13.d|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 <span id='citeF-26'></span>[[#cite-26|[26]]].

​

After discretization in space (Galerkin weighted residual method), and assuming <math display="inline">\Gamma _N</math> to be a free surface (i.e.&nbsp;<math display="inline">\bar \boldsymbol{\sigma }_n=\boldsymbol{0}</math>), the split discrete equations read:

​

<span id="eq-15"></span>

<span id="eq-16"></span>

<span id="eq-17"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>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 </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (15)

|-

| style="text-align: center;" | <math> 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 </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (16)

|-

| style="text-align: center;" | <math> 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} </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (17)

|}

|}

​

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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\boldsymbol{M}_{\rho }^{ab} = \int _{\Omega } \boldsymbol{N}^a \rho \boldsymbol{N}^b \;d\Omega </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (18)

|-

| style="text-align: center;" | <math>  \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 </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (19)

|-

| style="text-align: center;" | <math>  \boldsymbol{B}^{ab} = \int _{\Omega } (\boldsymbol{\nabla }\cdot \boldsymbol{N}^a)\; N^b \;d\Omega </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (20)

|-

| style="text-align: center;" | <math>  \boldsymbol{\boldsymbol{L}}_{(1/\rho )}^{ab} =\int _{\Omega } \boldsymbol{\nabla }N^a \frac{1}{\rho } \boldsymbol{\nabla }N^b \;d\Omega </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (21)

|-

| style="text-align: center;" | <math>  \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  </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (22)

|}

|}

​

where <math display="inline">\boldsymbol{N}</math> and <math display="inline">N</math> are the standard linear FE shape functions for velocity and pressure, and superscripts <math display="inline">a, b</math> 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.&nbsp;see discussion about numerical integration of discontinuous and singular functions in <span id='citeF-27'></span>[[#cite-27|[27]]]) or enrich the shape function space, as e.g.&nbsp;in <span id='citeF-28'></span><span id='citeF-29'></span>[[#cite-28|[28,29]]] and in the eXtended Finite Element Method (XFEM) <span id='citeF-30'></span><span id='citeF-31'></span>[[#cite-30|[30,31]]]. The nodal interface description for immiscible fluids in PFEM allows to evaluate exactly these integrals.

​

Conditions ([[#eq-12.c|12.c]]) and ([[#eq-12.d|12.d]]) are naturally included in <math display="inline">\boldsymbol{F}</math>, while pressure conditions ([[#eq-13.c|13.c]]) on <math display="inline">\Gamma _N</math> and ([[#eq-13.d|13.d]]) on <math display="inline">\Gamma _{int}</math> need to be weakly imposed in order to overcome the singularity of <math display="inline">\boldsymbol{\boldsymbol{L}}_{(1/\rho )}</math>:

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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  </math>

|}

|}

​

Both integrals are discretized in space as <math display="inline">\boldsymbol{M}_c (\boldsymbol{P}^{n+1}-\boldsymbol{P}^n)</math>, with <math display="inline">\boldsymbol{M}_c^{ab}=\int _\Gamma N^a N^b \, d\Gamma </math> the pressure mass matrix on the contours <math display="inline">\Gamma _N</math> and <math display="inline">\Gamma _{int}</math>, and incorporated into the pressure Poisson equation in the following way <span id='citeF-32'></span>[[#cite-32|[32]]]:

​

<span id="eq-23"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\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  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (23)

|}

​

The penalty parameter <math display="inline">\lambda </math> is defined as <math display="inline">\lambda =\alpha \;\frac{\Delta t}{h}</math>. It has to be sufficiently large so that the system matrix becomes invertible. <math display="inline">\alpha </math> is a scalar factor <math display="inline">\alpha \sim{\cal O}(1)</math> such that for <math display="inline">\alpha \rightarrow 0</math>, the pressure equation is satisfied but the discrete system may continue singular, and <math display="inline">\alpha \rightarrow \infty </math> 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 <math display="inline">{\cal C}^{-1}</math> discontinuity of the hydrostatic pressure gradient. Details are explained in <span id='citeF-2'></span>[[#cite-2|[2]]] and <span id='citeF-33'></span>[[#cite-33|[33]]].

​

==6 SURFACE TENSION==

​

The surface tension force at the interface, <math display="inline">\boldsymbol{f}_{st}=\int _{\Gamma _{int}} \gamma \kappa \boldsymbol{n}\cdot \boldsymbol{w}\,d\Gamma </math>, is naturally incorporated in the weak form of the finite element method. If it is discretized explicitly, i.e.&nbsp;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 <span id='citeF-34'></span>[[#cite-34|[34]]]:

​

<span id="eq-24"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\Delta t_{st}<\sqrt{\frac{\langle \rho \rangle \, h^3}{2\pi \gamma }}  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (24)

|}

​

where <math display="inline">\langle \rho \rangle =\frac{1}{2}(\rho _1+\rho _2)</math>, <math display="inline">h</math> is the mesh size and <math display="inline">\Delta t_{st}</math> the capillary time step. With this restriction the propagation of capillary waves is resolved and their unstable amplification avoided. Unfortunately, Eq.&nbsp;([[#eq-24|24]]) can be rather limiting for fine meshes and large surface tension coefficients. An implicit <span id='citeF-35'></span>[[#cite-35|[35]]] or semi-implicit <span id='citeF-36'></span>[[#cite-36|[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 <span id='citeF-34'></span>[[#cite-34|[34]]].

​

===6.1 Curvature calculation===

​

There are several ways to calculate the curvature <math display="inline">\kappa </math> 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 <math display="inline">\kappa \boldsymbol{n}</math> required for <math display="inline">\boldsymbol{f}_{st}</math> is calculated as:

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\boldsymbol{n}= \frac{\boldsymbol{R}}{|\boldsymbol{R}|}, \qquad \kappa \boldsymbol{n}= \frac{\boldsymbol{R}}{|\boldsymbol{R}|^2} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (25)

|}

​

<div id='img-5'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-curvature_computation2.png|300px|Calculation of the osculating circle at node x.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 5:''' Calculation of the osculating circle at node <math>\boldsymbol{x}</math>.

|}

​

===6.2 Static bubble===

​

We verify our surface tension algorithm with the static bubble test case <span id='citeF-37'></span><span id='citeF-26'></span>[[#cite-37|[37,26]]]. It consists in a circular fluid bubble into another viscous fluid at rest (see Figure [[#img-6|6]]), where gravitational or other external forces are neglected. According to the Laplace-Young law, the pressure jump will be <math display="inline">p_2-p_1=\gamma \kappa =\gamma /R</math>, where <math display="inline">p_2</math> is the bubble internal pressure, <math display="inline">p_1</math> the outer pressure and <math display="inline">R</math> 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. <div id='img-6'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-staticbubble_iniconfig.png|180px|Static bubble initial configuration.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 6:''' Static bubble initial configuration.

|}

​

We have simulated the equilibrium state of a circular bubble of a radius <math display="inline">R=0.25</math> (constant curvature <math display="inline">\kappa=1/R=4</math>) in a static fluid with <math display="inline">\rho _1 =\rho _2 =1</math>, <math display="inline">\mu _1 =\mu _2 =1</math>, <math display="inline">g=0</math>, and <math display="inline">\gamma=1</math>. The simulations have been run 100 time steps with <math display="inline">\Delta t</math> 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 [[#img-7|7]]) that may deform the bubble shape, produce a significant mass loss <span id='citeF-38'></span>[[#cite-38|[38]]] and spoil the solution.

​

<div id='img-7'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-staticbubble_parasiticvelocity2.png|186px|Spurious velocities in static bubble for continuous pressure.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 7:''' Spurious velocities in static bubble for continuous pressure.

|}

​

Figure [[#img-8|8]] shows the pressure solution for continuous approximation and different mesh sizes. The exact jump is <math display="inline">\Delta p = \gamma /R = 4</math>. 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 [[#img-9|9]], the solution is already excellent for coarse meshes. The velocity error (measured in the norm <math display="inline">||u||_{\infty }=\max _i |u_i|</math>) for both approximations is shown in Table [[#table-1|1]]. The discontinuous pressure produces velocity solutions three orders of magnitude better than the continuous one.

​

<div id='img-8a'></div>

<div id='img-8b'></div>

<div id='img-8'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-staticbubble_continuouspress_t1_comparison.png|294px|(a)]]

|[[Image:draft_Samper_466629376-staticbubble_contpressure.png|240px|(b)]]

|- style="text-align: center; font-size: 75%;"

| (a) 

| (b) 

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''Figure 8:''' Continuous pressure approximation: (a) profile at final time for different <math>h</math>, and (b) pressure field for <math>h=1/20</math>.

|}

​

<div id='img-9a'></div>

<div id='img-9b'></div>

<div id='img-9'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-staticbubble_discontinuouspress_t1_comparison.png|294px|(a)]]

|[[Image:draft_Samper_466629376-staticbubble_dispressure.png|240px|(b)]]

|- style="text-align: center; font-size: 75%;"

| (a) 

| (b) 

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''Figure 9:''' Discontinuous pressure approximation: (a) profile at final time for different <math>h</math>, and (b) pressure field for <math>h=1/20</math>.

|}

​

​

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"

|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1  Velocity error <math>||u||_{\infty }</math> at final time for continuous and discontinuous pressure approximations.

|- style="border-top: 2px solid;"

| <math display="inline">h</math> 

| Continuous 

| Discontinuous 

|- style="border-top: 2px solid;"

|  1/20 

| <math>4.4\times{10}^{-2}</math>

| <math>2.8\times{10}^{-5}</math>

|-

| 1/40 

| <math>1.7\times{10}^{-2}</math>

| <math>1.3\times{10}^{-5}</math>

|- style="border-bottom: 2px solid;"

| 1/80 

| <math>7.2\times{10}^{-3}</math>

| <math>8.9\times{10}^{-6}</math>

​

|}

​

Finally, Figure [[#img-10|10]] shows the influence of the parameter <math display="inline">\lambda </math> on the pressure solution at the interface. For the minimum value <math display="inline">\alpha _{min}</math> that makes the system Eq.&nbsp;([[#eq-23|23]]) invertible, the pressure profile is flat at the interface, i.e.&nbsp;the jump is represented exactly and incompressibility is satisfied. The larger the value of <math display="inline">\alpha </math>, the more strongly the pressure continuity condition at the interface is imposed.

​

<div id='img-10'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-staticbubble_discontinuouspress_t1_tauc_h005.png|276px|Pressure profile for discontinuous approximation and h=1/20. Influence of α value on the pressure jump at the interface: α<sub>min</sub>, 10α<sub>min</sub>, 20α<sub>min</sub>.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 10:''' Pressure profile for discontinuous approximation and <math>h=1/20</math>. Influence of <math>\alpha </math> value on the pressure jump at the interface: <math>\alpha _{min}</math>, <math>10\alpha _{min}</math>, <math>20\alpha _{min}</math>.

|}

​

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 <span id='citeF-39'></span>[[#cite-39|[39]]] developed a well-known graphical correlation for single gas bubbles rising in an infinite liquid using three dimensionless numbers: the Reynolds number (<math display="inline">Re</math>), the Eötvös number (<math display="inline">Eo</math>), and the Morton number (<math display="inline">Mo</math>).

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>Re = \displaystyle \frac{\rho \, L \, U}{\mu } </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (26)

|-

| style="text-align: center;" | <math> Eo = \displaystyle \frac{\rho \, g \, L^2}{\gamma } </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (27)

|-

| style="text-align: center;" | <math> Mo = \displaystyle \frac{g \, \mu ^4}{\rho \, \gamma ^3} </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (28)

|}

|}

​

where <math display="inline">L</math> is an equivalent diameter of the bubble, and <math display="inline">U</math> the rising speed of the bubble. <span id='citeF-39'></span>[[#cite-39|[39]]] classifies the bubble shapes into three regimes: spherical, ellipsoidal, and spherical cap. Bubbles with low <math display="inline">Re</math> or low <math display="inline">Eo</math> are spherical. For higher <math display="inline">Re</math>, bubbles have an ellipsoidal shape at intermediate <math display="inline">Eo</math> and spherical cap shape at high <math display="inline">Eo</math>. A more detailed regime diagram was given by Clift et al.&nbsp;<span id='citeF-40'></span>[[#cite-40|[40]]], which included wobbling bubbles for <math display="inline">Re \sim 10^3</math> in the ellipsoidal regime, and the spherical cap regime was subdivided into spherical cap, skirted, and dimpled ellipsoidal cap for high, intermediate, and low <math display="inline">Re</math> numbers, respectively. These diagrams were further developed by Bhaga et al.&nbsp;<span id='citeF-41'></span>[[#cite-41|[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 <math display="inline">Re</math>.

​

Numerous experimental studies have been performed to understand the flow dynamics of a single rising bubble, e.g.&nbsp;by <span id='citeF-42'></span><span id='citeF-39'></span><span id='citeF-43'></span><span id='citeF-40'></span><span id='citeF-41'></span><span id='citeF-44'></span><span id='citeF-45'></span><span id='citeF-46'></span>[[#cite-42|[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  <span id='citeF-47'></span>[[#cite-47|[47]]], <span id='citeF-48'></span>[[#cite-48|[48]]], and <span id='citeF-49'></span>[[#cite-49|[49]]] used the front-tracking method (together with finite differences), also <span id='citeF-50'></span>[[#cite-50|[50]]] (finite differences) and <span id='citeF-51'></span>[[#cite-51|[51]]] (finite volume method); level set is used by <span id='citeF-27'></span><span id='citeF-26'></span><span id='citeF-1'></span>[[#cite-27|[27,26,1]]] (finite elements) and <span id='citeF-52'></span>[[#cite-52|[52]]] (finite volumes); Volume-of-Fluid is used in <span id='citeF-53'></span>[[#cite-53|[53]]] and <span id='citeF-54'></span>[[#cite-54|[54]]]; <span id='citeF-55'></span>[[#cite-55|[55]]] use a hybrid approach between VOF and level set (finite volumes); and interface fitting method in <span id='citeF-56'></span>[[#cite-56|[56]]] and <span id='citeF-45'></span>[[#cite-45|[45]]]. The Lattice-Boltzmann method is used e.g.&nbsp;in <span id='citeF-57'></span>[[#cite-57|[57]]] and <span id='citeF-58'></span>[[#cite-58|[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.&nbsp;<span id='citeF-1'></span>[[#cite-1|[1]]]. In Section [[#7.1 Comparison with previous numerical experiments|7.1]] we compare the PFEM results with these test solutions, and in Section [[#7.2 Bubble breakup and coalescence|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===

​

<div id='img-11'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-confini.png|180px|Initial configuration.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 11:''' Initial configuration.

|}

​

The problem consists in a bubble rising in a liquid column as illustrated in Figure [[#img-11|11]]. Two tests have been proposed in <span id='citeF-1'></span>[[#cite-1|[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 [[#table-2|2]].

​

​

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"

|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Physical parameters.

|- style="border-top: 2px solid;"

|  Test case 

| <math>\rho _1</math>

| <math>\rho _2</math>

| <math>\mu _1</math>

| <math>\mu _2</math>

| g 

| <math>\gamma </math>

| <math>Re</math>

| <math>Eo</math>

| <math>\rho _1/\rho _2</math>

| <math>\mu _1/\mu _2</math>

|- style="border-top: 2px solid;"

|  1 

| 1000 

| 100 

| 10 

| 1 

| 0.98 

| 24.5 

| 35 

| 10 

| 10 

| 10

|- style="border-bottom: 2px solid;"

| 2 

| 1000 

| 1 

| 10 

| 0.1 

| 0.98 

| 1.96 

| 35 

| 125 

| 1000 

| 100

​

|}

​

The reference solutions presented in <span id='citeF-1'></span>[[#cite-1|[1]]] have been run with three different numerical approaches: the TP2D <span id='citeF-59'></span>[[#cite-59|[59]]], the FreeLIFE <span id='citeF-60'></span>[[#cite-60|[60]]], and the MooNMD <span id='citeF-61'></span>[[#cite-61|[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 <math display="inline">t=3</math> s,

* circularity <math display="inline">\displaystyle c=\frac{2\sqrt{\pi \mbox{Area}}}{\mbox{Perimeter}}</math>,

* center of mass <math display="inline">\displaystyle \boldsymbol{X}_c=\int _{\Omega _2}\boldsymbol{x}\,dx/\int _{\Omega _2}1\,dx</math>, and

* rise velocity <math display="inline">\displaystyle \boldsymbol{U}_c=\int _{\Omega _2}u\,dx/\int _{\Omega _2}1\,dx</math>.

​

The computations have been performed on unstructured meshes with element size <math display="inline">h=1/40</math> in the bulk of the fluids and wall regions, and the mesh at the interface has been refined to <math display="inline">h=1/80</math>, <math display="inline">1/160</math>, <math display="inline">1/320</math> and <math display="inline">1/640</math> in order to analyze the convergence in <math display="inline">h</math> of the solution. With a refinement based on the distance to the interface (see Figure [[#img-12|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.

​

<div id='img-12'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-bubble_meshref1_bw.png|180px|]]

|[[Image:draft_Samper_466629376-bubble_meshref2_bw.png|240px|Mesh is refined close to the interface with the help of a distance function.]]

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''Figure 12:''' Mesh is refined close to the interface with the help of a distance function.

|}

​

<div id='img-13'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|

{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"

|-

| [[Image:draft_Samper_466629376-test1_bubble_mat_bw-01.png|300px|figures/test1_bubble_mat_bw-01.png]]

| [[Image:draft_Samper_466629376-test1_bubble_mat_bw-02.png|300px|figures/test1_bubble_mat_bw-02.png]]

| [[Image:draft_Samper_466629376-test1_bubble_mat_bw-03.png|300px|figures/test1_bubble_mat_bw-03.png]]

| [[Image:draft_Samper_466629376-test1_bubble_mat_bw-04.png|300px|figures/test1_bubble_mat_bw-04.png]]

|-

| (a) <math display="inline">t=0</math> s 

| (b) <math display="inline">t=0.5</math> s 

| (c) <math display="inline">t=1.0</math> s 

| (d) <math display="inline">t=1.5</math> s 

|-

| [[Image:draft_Samper_466629376-test1_bubble_mat_bw-05.png|300px|figures/test1_bubble_mat_bw-05.png]]

| [[Image:draft_Samper_466629376-test1_bubble_mat_bw-06.png|300px|figures/test1_bubble_mat_bw-06.png]]

| [[Image:draft_Samper_466629376-test1_bubble_mat_bw-07.png|300px|figures/test1_bubble_mat_bw-07.png]]

| 

|-

| (e) <math display="inline">t=2.0</math> s 

| (f) <math display="inline">t=2.5</math> s 

| (g) <math display="inline">t=3.0</math> s 

|}

​

|-

​

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 13:''' Test 1. Bubble evolution for mesh size <math>h=1/320</math>.

|}

​

For test 1, Figure [[#img-13|13]] shows the evolution of the rising bubble. At final time <math display="inline">t=3</math> s we compare the PFEM bubble shapes for the meshes <math display="inline">h=1/40, 1/80, 1/160</math> and <math display="inline">1/320</math>, and observe that they converge to the shape of the finest mesh (Figure [[#img-14|14]]), which is in good agreement with the TP2D solution reported in <span id='citeF-1'></span>[[#cite-1|[1]]] (Figure [[#img-15|15]]a). The plots of the bubble circularity (Figure [[#img-15|15]]b) and rise velocity (Figure [[#img-15|15]]d) show that our bubble is slightly oscillating, but the evolution of the center of mass (Figure [[#img-15|15]]c) 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, <math display="inline">\Delta V=\frac{V_f-V_0}{V_0}</math>, is of order <math display="inline">{\cal O}(10^{-4})</math> (Table [[#table-3|3]]).

​

<div id='img-14'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_466629376-test1_bubbleshape-h.png|336px|Test 1. PFEM bubble shape for different mesh sizes at t=3 s.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 14:''' Test 1. PFEM bubble shape for different mesh sizes at <math>t=3</math> s.

|}

​

<div id='img-15a'></div>