A stabilised ﬁnite element framework for viscoelastic multiphase ﬂows using a conservative level-set method

The development of stable numerical schemes for viscoelastic multiphase ﬂows has implications for many areas of engineering applications. The principal original contribution of this paper is the implementation of a conservative level-set method to deﬁne implicitly the interface between ﬂuid phases, fully integrated into the mathematical framework of viscoelastic ﬂow. The governing equations are discretized using the ﬁnite element method and stabilisation of the constitutive equation is achieved using either the discontinuous Galerkin (DG) or streamline upwinding (SU) method. The discrete elastic viscous stress splitting gradient (DEVSS-G) formulation is also employed in the Navier-Stokes equations to balance the hyperbolic characteristics of the polymeric stress tensor. The numerical scheme is validated with reference to several benchmark problems and excellent quantitative agreement with published data is found for Newtonian and viscoelastic ﬂuids, for both single and multiphase ﬂows. The motion of a gas bubble rising in a viscoelastic ﬂuid is studied in detail. The inﬂuence of polymer concentration, surface tension, ﬂuid elasticity and shear-thinning behaviour, on ﬂow features such as the development of ﬁlaments and cusps and the generation of negative wakes is explored.  2023 The Author(s). Published by Elsevier Inc. This is an open access article under the


Introduction
Multiphase flows are ubiquitous in a variety of natural and industrial processes, from the refinement of oil to food science and the flow inside a nuclear reactor.Furthermore, considering one of the phases as having a viscoelastic relationship between stress and strain, rather than a traditional Newtonian relationship, is important in chemical manufacturing and polymer processing.Understanding how these complex flows evolve over time is of great importance.Using numerical schemes, like the one presented in this paper, rather than performing physical experiments reduces cost and gives us more information over what parametric changes cause the phenomenological activity we observe.
The main object that requires rigorous definition when studying multiphase flows is the interface between the fluid phases.The methods used can be broadly split into two categories: interface capturing (e.g.volume of fluid (VOF) [1] and level-set [2]) and interface tracking (e.g.front tracking [3] and moving mesh method [4]).Interface tracking methods explicitly construct a set of points that define the interface which are moved according to the local velocity.Interface capturing methods only implicitly define the interface, for example, as the contour of a greater surface or the amount of marker in a cell.In this paper we opt for an implicit representation of the interface since viscoelastic fluids can have very unpredictable dynamics, so the ability to track the evolution of very general interfaces is desirable.The level-set method uses the 0-contour of a signed distance function to define the interface, which is then advected in a velocity field and reinitialised to maintain its original properties.Initially introduced by Osher and Sethian [5], it was eventually applied to incompressible two-phase flow by Sussman et al. [2], with applications ranging from topology optimization [6] to mathematical biology [7].The presence of steep polymeric stress boundary layers on the interface can contribute to the loss of mass from one phase over time.To remedy the lack of mass conservation of the traditional level-set method, we utilise the conservative level-set method of Olsson et al. [8].To the best of our knowledge this has not been applied to viscoelastic flow before.The conservative level-set method uses a smeared Heaviside function rather than a signed distance function in order to represent the interface.The reinitialisation scheme now involves taking an L 2 projection of the normal to the interface and compressing the surface in that direction.This modification of the traditional level-set method has resulted in many new developments, with Kees et al. [9] adapting it for variable-order approximations / unstructured meshes and Zahedi et al. [10] using it for contact line dynamics.More recent developments of the method come from Lin et al. [11] who used the conservative level-set method for modelling multiphase thermo-fluid flow and Jettestuen et al. [12] who modelled capillary-controlled displacements in porous media.
When accounting for viscoelastic dynamics, the problem gets more challenging on multiple levels.Firstly, we have to solve a separate constitutive equation that relates polymeric stress and deformation rate.This equation is hyperbolic, unlike the conservation of mass and momentum equations which form an elliptic system.This mandates the use of numerical stabilisation techniques in order to obtain converged approximations to our problem.One technique is a version of the popular elastic viscous split stress (EVSS) [13,14] type methods, namely the discrete EVSS gradient (DEVSS-G) [15] method.Here, an L 2 projection of the velocity gradient tensor is used to approximate velocity gradient terms in the momentum and constitutive equations.This has the effect of creating artificial ellipticity, while maintaining consistency with the original equations.The discontinuous Galerkin (DG) [16] finite-element discretisation uses discontinuous finite elements for polymeric stress to better approximate and resolve the steep stress gradients at boundary layers.Through the choice of the numerical flux function, we can also ensure a fully upwind approximation in the process.The last stabilisation technique is the streamline upwinding (SU) [17] method, which modifies the test functions in the weak formulation, introducing artificial diffusion into the constitutive equation in the process.Numerical methods for solving viscoelastic flow problems are well-established.The review by Baaijens [18] describes many of the popular stabilisation methods like DG, SU and DEVSS-G in detail.
There have been a number of contributions to the development of numerical methods for viscoelastic multiphase flows in recent years.A short literature review is provided here of some of the most important papers.Figueiredo et al. [19] investigated the use of the VOF method to model the famous Weissenberg effect, while Izbassarov et al. [20] used the front tracking method to investigate viscoelastic emulsions.Harvie et al. [21] considered a VOF approach for modelling a Newtonian droplet passing through a microfluidic contraction and Chung et al. [22] considered a similar 5:1:5 contraction/expansion microchannel using a finite-element front-tracking approach.Vahabi et al. [23] and Zainali et al. [24] investigated the viscoelastic rising bubble problem by using different versions of the smoothed particle hydrodynamics (SPH) [25] method.Pillapakkum and Singh [26] developed a level-set method for modelling a Newtonian droplet in a viscoelastic fluid under simple shear, with a brief investigation into the rising bubble problem as well.Pillapakkum et al. [27] expanded on this work some years later and constructed a fully 3D model for a rising bubble in a viscoelastic fluid, also based on the level-set method.Venkatesan et al. [28] considered an arbitrary Lagrangian-Eulerian (ALE) approach for a Newtonian bubble rising in a viscoelastic bulk fluid and a viscoelastic bubble rising in a Newtonian bulk fluid.Lind and Phillips [29] developed a boundary element method approach to model a rising Newtonian bubble in a viscoelastic fluid, which included a detailed parametric study.Recently a lot of attention has focused on modelling the jump discontinuity in the rise velocity of a rising gas bubble in a viscoelastic fluid, demonstrated experimentally in Pilz and Brenn [30].This phenomenon is a sharp jump in the rise velocity as a function of bubble volume, where there exists a so-called critical volume at which the rise velocity increases drastically.Yuan et al. [31] and Niethammer et al. [32] used an interface capturing approach through the VOF method to construct fully 3D simulations for modelling the velocity jump discontinuity.Both papers accounted for such large Weissenberg numbers that their numerical schemes necessitated the use of the log-conformation formulation [33] of the constitutive equation.Finally, Fraggedakis et al. [34] adopted a theoretical approach instead and performed a steady state analysis with the aid of pseudo arc length computations to model the jump discontinuity.
To summarise the novelty of our paper, the conservative level-set method is applied to viscoelastic flow in a 2D finite element formulation with stabilisation and non-dimensionalisation procedures included.This has only previously been done for the non-conservative level-set method, see Pillapakkam et al. [26] as an example.The conservative level-set method is more challenging to implement as the smeared Heaviside function is harder to approximate for sharp interfaces and requires frequent L 2 projections for normal calculations.The reward is a numerical scheme that accurately calculates the mass of a viscoelastic rising bubble while it undergoes strong polymeric stress forces near the interface.The viscoelastic constitutive equation is stabilised using a variety of intricate techniques in order to accurately model flows where the extension rate is large.Papers from Vahabi et al. [23] and Zainali et al. [24] focus only on small polymeric concentrations and thus do not require any stabilisation, so we expand on this work through the aforementioned methods.We validate our results using benchmarks from the literature, obtaining very satisfactory results for both the multiphase Newtonian and singlephase viscoelastic cases.Lastly, the case of a Newtonian bubble rising in a viscoelastic fluid is explored.We use additional viscoelastic models when compared with the literature and model much larger polymer concentrations, which allows us to analyse phenomena in detail.A fully object-oriented, parallelised and open-source code based on the FEniCS [35] finite element library is available in the lead author's Github [36] repository for reproduction of all results found in the paper.
The paper is organised as follows.In Section 2 the mathematical model for describing viscoelastic multiphase flows is presented.In Section 3 the finite element discretization of the governing equations is described paying particular attention to the use of stabilisation techniques.The results of numerical experiments on several benchmark problems are presented and discussed in Section 4. Finally, in Section 5 some concluding remarks and areas of future study are summarised.

Multiphase governing equations
The multiphase governing equations are comprised firstly of the following mass and momentum conservation laws: where u is the velocity vector, ρ is the density and F is a body force.In these equations σ is the Cauchy stress tensor defined by: where p is the pressure, η s is the solvent viscosity (η s = μ for Newtonian fluids), D (u) = 1 2 ∇u + ∇u T is the rate-of-strain tensor and τ p is the polymeric stress tensor.For Newtonian fluids, τ p = 0.For viscoelastic fluids the evolution of τ p is dictated by a constitutive equation, detailed more in Section 2.2.In the momentum equation the body force F comprises of a gravitational force F g and a surface tension force F : The treatment of the surface tension force follows the continuum surface force approach of Brackbill et al. [37] where the force due to surface tension is transformed into a volume force.Here σ is the surface tension coefficient, φ is the conservative level-set function (see Section 2.3) and κ is the curvature of the interface.We refer to Section 2.4 for more information on the continuum surface force approach.The material parameters μ and ρ are now phase-dependent: where μ i and ρ i are the values of the dynamic viscosity and density in phase i. Eqs.(1) are solved in a time dependent twodimensional domain (t) = 1 (t) ∪ 2 (t) where i (t) is the domain occupied by phase i at time t ∈ (0, T ].The two phases are distinguished by the interface between them (t), which itself is determined from the 0.5 contour of the conservative level-set function.A schematic figure of this set-up is shown in Fig. 1.

Viscoelastic flow
For a viscoelastic fluid the evolution of the polymeric stress tensor τ p can be described by a multitude of different constitutive laws.Cracco [38] provides a general form of the constitutive equation, from which multiple different relationships can be obtained by setting certain parameters to 0 or 1.In this paper we consider the following relationship: where τ p denotes the upper-convected derivative of the polymeric stress τ p defined by: For a viscoelastic fluid, the total viscosity η 0 is defined as the sum of two viscosities: the solvent viscosity η s and the polymeric viscosity η p .Both of these material parameters can be obtained by considering the amount (by weight) of polymer dissolved in the Newtonian solvent, denoted by c.The relaxation time λ 1 is the characteristic time taken for the shear stress in a simple shear flow to relax under constant strain conditions.The parameter α is the Giesekus mobility factor and controls the shear-thinning behaviour of the fluid.When using the Oldroyd-B model [39], we shall use α = 0. Oldroyd-B is largely unsuitable for extensional flows due to the presence of a singularity in the extensional viscosity η E at a finite elongation rate ˙ .In contrast the Giesekus model predicts a bounded extensional viscosity.To illustrate this, consider a steady uniaxial elongational flow field of the form u = ˙ x, − ˙ 2 y, − ˙ 2 z .The extensional viscosity is given by the expression: where τ xx and τ yy are the normal components of the polymeric stress tensor τ p .For the given uniaxial extensional flow field the constitutive equation ( 5) can be solved for each value of ˙ to obtain the corresponding value of η E ( ˙ ).The extensional viscosity is plotted in Fig. 2 for the Oldroyd-B and Giesekus models, the latter for different values of the mobility factor α.

Level-set method
The level-set function used for the conservative level-set method is generated from a transformation of the signed distance function where φ (x) this defines the interface implicitly.We perform the following transformation to obtain the level-set function used in the conservative level-set method: This transforms the signed distance function into a smeared Heaviside function, with the interface located at φ(x) = 0.5.
The parameter is a measure of interface thickness and determines the width of the transition region between the 0 level-set and 1 level-set.To capture the interface we solve the advection equation for φ with a given velocity field u.For incompressible fluids we have ∇ • u = 0, so it can be shown the advection equation is equivalent to its conservative form: Hence the advection problem consists of solving Eq. ( 9) for φ with initial condition φ(x, 0) = φ 0 .During advection of the level-set function, the desirable stability properties of a smeared Heaviside function are lost.Therefore an intermediate reinitialisation step that recaptures these properties must be performed, while also not distorting the location of the interface.Upon solving Eq. ( 9) for an intermediate φ we calculate the L 2 -projection of the interface normal n and reinitialise by solving the following equation: until a steady state is reached under the artificial time-step τ .On the continuous level, the interface normal is given by n = ∇φ |∇φ| .Since Eq. ( 10) is written in conservative form we know that φ is also conserved as it is reinitialised, which means the area bounded by the φ (x) = 0.5 contour is conserved under advection and reinitialisation.This proves that the method conserves the mass of the level-set function.The area bounded by the φ = 0.5 contour is also conserved, the proof of which can be found in Olsson et al. [40].

Boundary conditions
Eqs. (1) are solved subject to prescribed initial conditions: u(x, 0) = u 0 (x) (11) and boundary conditions, which can be of Dirichlet: (12) or Neumann type: The terms t and n refer to the tangent and normal vectors to the interface respectively.Along the interface (t) we enforce the following conditions: (16) where we drop x and t for brevity.The first of the interface conditions (Eq. ( 15)) ensures that the velocity field is continuous across (t) while the second (Eq.( 16)) states that the jump in the normal stress is balanced by the force due to surface tension.The second condition is satisfied by the continuum surface force approach of Brackbill et al. [37], in which surface tension is interpreted as a continuous volume force in the neighbourhood of the interface (F in Eq. ( 3)).This is accomplished by expressing Eq. ( 16) as a volume force: (17) where δ( (t)) is the Dirac delta function localised to the interface and κ is the curvature given by κ = ∇ • n .Then using the fact that δ( (t)) = |∇φ| we derive the form of F given in Eq. (3).We refer to Brackbill et al. [37] and Chang et al. [41] for a full derivation.The approach is also implemented in the papers cited in Sections 4.1 and 4.3.These include Hysing et al. [42], Vahabi et al. [23] and Zainali et al. [24].

Numerical method
The finite element method (FEM) is used to discretise the governing equations presented in Section 2.1.The FEM is based on the weak formulation of the governing equations.The computational domain is decomposed into finite elements on which a piecewise polynomial approximation is sought.

Function spaces and weak formulations
In order to define the weak formulation of the problem, we first define suitable function spaces that guarantee the wellposedness of the integrals that appear in the weak forms.For the level-set advection and reinitialisation we choose the trial space φ ∈ H 1 ( ) and test space ψ ∈ H 1 0 ( ).For velocity we choose the trial space u ∈ H 1 D ( ) , where the subscript D denotes satisfaction of given Dirichlet boundary conditions, and corresponding test space v ∈ H 1 0 ( ) Here d represents the dimension of the problem.Since pressure only appears as a gradient term in the Navier-Stokes equations it is only defined up to an arbitrary constant.To guarantee the uniqueness of pressure we choose the following trial and test spaces: Lastly, the polymeric stress is a second order symmetric tensor so its trial and test functions τ p and S must reflect this.
Define the following inner products for scalar, vector and tensor functions in L 2 respectively: The weak formulation of Eqs. ( 1) and ( 5) is: For t ∈ (0, T ], we find (u, p, τ p ) ∈ d×d sym with u(x, 0) = u 0 and τ p (x, 0) = τ p,0 .We group the body force terms under a single operator f : and the trilinear operator a and bilinear operators b and c are defined as follows: The weak formulations for Eqs. ( 9), (10) and the normal projection are as follows: For t ∈ (0, T ], we find φ(t) ∈ H 1 ( ) such that: where φ is the solution to Eq. ( 25) and δ is some small quantity which is added to the denominator to avoid division by 0. Finally, we find the reinitialised level-set function φ(τ ) ∈ H 1 ( ) from the following: ∀ψ ∈ H 1 0 ( ).The initial condition for Eqs. ( 26) and ( 27) is the solution of Eq. ( 25).Eq. ( 27) is solved until steady state is reached (see Section 3.5).

Galerkin method
Following the construction of the weak forms in Section 3.1 we formally state the finite-element approximation to the problem.First the computational domain is partitioned into a finite set of triangular elements T with the following properties: where ∅ is the empty set and i, j ∈ N. The non-boundary vertices of each cell will have a node with position vector x i for 1 ≤ i ≤ N where N is the number of nodes.Continuous piecewise linear basis functions ξ i with the property ξ i (x j ) = δ ij are also defined.Consider a finite-dimensional subspace U h of one of our infinite spaces (say H 1 D ( )) from earlier.The basis functions ξ i , 1 ≤ i ≤ N, span this subspace and as a result u h ∈ U h can be represented as: (29) where the coefficients U i are the velocity degrees of freedom.For the Navier-Stokes equations, the finite element approximation to the weak formulation is to find be replaced by basis functions, which then allows the integrals to be evaluated and matrices to be constructed.This approximation is repeated for all the weak forms in Section 3.1.

Choice of finite elements
In order to guarantee a stable solution to the FE approximation of the Navier-Stokes equations, we must choose the pressure and velocity spaces in accordance with the Ladyzhenskaya-Babuška-Brezzi (LBB) condition.
The finite dimensional subspaces U h and P h are spaces of continuous piecewise linear polynomials, or in set theoretic notation: where P k refers to the space of polynomials of degree less than or equal to k. Pairs of finite elements with polynomial degrees P k+1 − P k satisfy the LBB condition and are known as Taylor-Hood elements.The most common choice, which is used in this paper, is the P 2 − P 1 velocity / pressure coupling (Fig. 3).The polymeric stress is approximated using P 2 elements.In Sections 3.4.1 and 3.4.2we will discuss the stabilisation techniques that are required in flows with high Weissenberg numbers.We found P 1 elements were insufficient for resolving the 0.5-level-set so P 2 elements were used for all level-set related equations.The choice of mesh resolution will be discussed in the results section but structured meshes are used throughout the paper.

Viscoelastic stabilisation techniques
In order to stabilise the viscoelastic constitutive equation, we introduce three different methods that will be used throughout this paper.For high Weissenberg numbers and large polymer concentrations the numerical solution to the constitutive equation grows exponentially and some form of stabilisation is necessary to obtain converged approximations.

Discontinuous Galerkin
The viscoelastic constitutive equation is a non-linear hyperbolic PDE which can give rise to steep stress gradients or thin boundary layers.For high Weissenberg number flows these regions of large polymeric stress become more prevalent in regions of high elastic deformation.One way to increase stability is to use a discontinuous Galerkin treatment of the boundary terms.This relaxes the continuity requirements described in Section 3.3 and allows the polymeric stress to be discontinuous across element boundaries.Integrate the convective term in the constitutive equation by parts over each element T i and sum over all elements in the triangulation T : The boundary integral is now expressed in terms of 'positive' and 'negative' contributions (Fig. 4), which refers to the two perspectives of the same cell boundary ∂ T i , since τ p is dual-valued there.The last term on the right hand side of Eq. ( 35) is reformulated as: where τ + p and τ − p are polymeric stresses evaluated on either side of their shared boundary.These will be distinct since the polymeric stress is discontinuous across elements.Of course it is true that some elements have unshared boundaries, namely those that have an edge which coincides with the Dirichlet boundary.We omit the unshared boundary term in Eq. ( 36) however as the boundary conditions we impose cause it to vanish anyway.In order to define the transport of τ p across boundaries, the term u • n ± τ ± p is approximated by a quantity known as the numerical flux.The numerical flux can be defined according to the characteristics of the conservation law under inspection, so if we expect information to be propagated downstream using information from upstream, we use a fully upwind approach.Many more choices of numerical flux are also possible.The numerical flux F is of the following form: where γ = 1 guarantees a fully upwind approach.Following the work of Pietro et al. [43], we use the upwind numerical flux that is formulated in terms of jumps and averages instead, which is equivalent to Eq. ( 37) for γ = 1: where {τ p } = We are now left with the following expression: We substitute Eq. ( 40) back into the weak form of the constitutive equation, then integrate by parts the first term on the right hand side of Eq. ( 40).This will retrieve the original advective term and a boundary term vanishes due to the incompressibility condition.Altogether this yields the following discontinuous Galerkin approximation to the constitutive equation:

Streamline upwinding
An alternative to DG is the streamline upwinding (SU) method [44].This method augments the test functions with an artificial diffusion term in order to stabilise the convective term: where the SU parameter γ will be specified in Section 4.2.The method was introduced and applied to the incompressible Navier-Stokes equations by Brooks and Hughes [17] and to the viscoelastic constitutive equation by Marchal and Crochet [44].Due to the presence of oscillatory stress fields near steep stress boundary layers Marchal and Crochet [44] applied the augmented test function to the convective term only.Even though this results in an inconsistent formulation and we are left with a residual, the increase in the reliability of the scheme is significant.The resulting weak form of the constitutive equation is as follows:

DEVSS-G
The last stabilisation method that is employed is a variant of the EVSS class of techniques that introduces artificial ellipticity to the momentum equation.A summary of this method is available in Baaijens [18] and has been considered for the case of viscoelastic multiphase flow by Chung et al. [22].The momentum equation includes forces due to Newtonian stress through the term η s ∇ • D(u) and forces due to elastic stress through the term ∇ • τ p .For highly elastic viscoelastic fluids, the first term is small when compared to the latter, leading to a reduced elliptic contribution in the momentum equation.One way to circumvent this is to consider the following approximation of the Cauchy stress tensor: where η θ = 1 − η s .In Eq. (44) we add the difference between two different approximations to the velocity gradient tensor, thereby increasing the ellipticity.In the continuous limit, the L 2 projection G calculated in Eq. ( 45) is equivalent to ∇u so the extra terms cancel.However, at the discrete level the projection is not exact and we have extra elliptic terms in the momentum equation, which help balance the hyperbolic characteristics created through the inclusion of polymeric stress.

Temporal discretisation
We use 2 nd order Crank-Nicolson and 1 st order implicit schemes for the advection and reinitialisation of the levelset function respectively.The conservation of mass and momentum equations are linked temporally through the use of a projection algorithm.The constitutive equation follows similarly to the level-set equations, using a 2 nd order Crank-Nicolson scheme.

Navier-Stokes equations
Due to the presence of the incompressibility condition and the fact that there is no equation of state for pressure, we must find a way of coupling the Navier-Stokes equations so as to produce a consistent numerical method.In this paper we use a variation on Chorin's classic projection method, the Incremental Pressure Correction Scheme (IPCS) of Goda [45].In this scheme we first neglect the incompressibility condition and solve the momentum equation to determine an intermediate velocity u * : The original method used a time-averaged velocity U = 1 2 (u * + u n ) in the viscous term but we found better convergence with the fully implicit choice of u * .
Then, solve a Poisson equation for the pressure correction, simultaneously ensuring ∇ • u n+1 = 0: Lastly, update the velocity field:

Constitutive equation
We use a 2 nd order Crank-Nicolson scheme for temporal discretisation: where T n is defined as: and the trilinear operators M and N are defined as follows:

Level-set equations
For the level-set advection we a 2 nd order Crank-Nicolson scheme as well: It is important to note again that the normal projection step in Eq. ( 26) comes immediately before the reinitialisation step and not simultaneous with it, meaning the same normal vector is used throughout the reinitialisation.A simple 1 st order implicit scheme was found to be suitable for advancing the reinitialisation equation through a pseudo time-step defined by where Eq. ( 54) is solved recursively until the following condition is satisfied: where δ is chosen to be some small quantity.

Non-dimensionalisation
Here we summarise the dimensionless form of the Navier-Stokes equations used in this paper.A similar nondimensionalisation procedure to that in Venkatesan et al. [28] is used with the following scales and dimensionless parameters: where Fr is the Froude number, Re is the Reynolds number, Wi is the Weissenberg number, Eo is the Eötvös number and β and η are viscosity ratios useful in the non-dimensionalisation.In our simulations Fr = 1 so U = √ gL where U is the velocity scale, L is the length scale and g is the acceleration due to gravity.The equations are scaled relative to phase 2. The total viscosity of the fluid is η 0,2 = η s,2 + η p,2 which is the sum of the solvent and polymeric viscosities respectively.The non-dimensional equations governing the flow of an incompressible viscoelastic fluid are (bar over dimensionless quantities has been removed): where F = ρ Fr 2 g − 1 Eo ∇φ.The dimensional density ρ remains in the dimensionless momentum equation but is now a nondimensional density ratio, scaled according to phase 2. Using this set of dimensionless equations, we construct a new weak from that includes the stabilisation techniques introduced in Section (3.4).The set of equations that follow are used for all results generated in Section 4. For t ∈ (0 d×d sym with u(x, 0) = u 0 , p(x, 0) = p 0 and τ p (x, 0) = τ p,0 .We also have the body force terms grouped under a single operator:

κ∇φ, v
The term denoted by Stab refers to the additional stabilising term that is unique to either the DG or SU method:

Table 2
The meshes used in Section 4.1 with their associated number of cells, minimum cell length and degrees of freedom (DOF) for each variable.The Mc mesh is not used in any computations but is provided for illustration purposes to show the structure of the mesh in Fig. 5. the limited available work from the literature, and work with much larger polymer concentrations too.We then perform a parametric study, exploring how the dynamics of the viscoelastic fluid and Newtonian bubble change with the relative importance of dimensionless quantities, namely the Eötvös and Weissenberg numbers, and the Giesekus mobility factor.

Benchmark: two-phase Newtonian flow
We validate our Newtonian multiphase model by making comparisons with the wealth of available literature for the 2D rising bubble benchmark case.In particular, the paper by Hysing et al. [42] reports on the findings of three different pieces of software.The 'TP2D' and 'FreeLIFE' codes employ an Eulerian level-set method while the 'MooNMD' code used an arbitrary Lagrangian-Eulerian (ALE) approach.We omit the 'TP2D' data as numerically it is effectively the same as 'FreeLIFE' but its quantitative behaviour, especially for Case 2, is very different.The computational domain is (t) = [0, 1] × [0, 2].The bubble's initial position is x i = (0.5, 0.5) and has radius r = 0.25.Slip boundary conditions are employed on the left and right walls and no-slip conditions on the top and bottom walls: where D and N correspond to the solid boundaries labelled in Fig. 5. Interface boundary conditions ( 16) are enforced by their inclusion in the weak formulation.The relevant material parameters can be found in Table 1.The interface thickness and reinitialisation pseudo time-step are chosen to be = 1.5 x and τ = 0.1 x respectively, where x is the minimum cell size.Four different bubble metrics: shape, centre of mass x c , rise velocity U r and circularity C are used for quantitative comparison purposes: We also compare the final bubble shape against the reference data.The relevant material parameters are given in Table 1.Three different meshes are used, each with increasing levels of refinement, detailed in Table 2.

Case 1
Fig. 6 compares the predictions obtained using the scheme described in this paper with 2 of the 3 sets of benchmark data provided by Hysing [42] for Case 1. Excellent agreement is observed with the benchmark results and convergence is attained with all of the bubble metrics.Quantitatively the bubble forms an ellipse at around T = 1, with the cap widening until around T = 2. Subsequently the bubble maintains this shape while continuing to rise for the rest of the simulation time.Due to the small Eötvös number, surface tension effects are small and the bubble undergoes only moderate interfacial deformation in response to buoyancy forces.

Case 2
The corresponding predictions of the bubble metrics for Case 2 (see Table 1) are presented in Fig. 7.The Reynolds and Eötvös numbers are both larger than in Case 1.This bubble lies within the ellipsoidal cap regime which is characterised by the bubble inverting in on itself, then surface tension forces acting to correct the bubble shape while thin trailing filaments are left behind in the wake.Bubble breakup can be exhibited but as for the 'MooNMD' code, this behaviour is not predicted.Similar to Case 1, the centre of mass and rise velocity metrics agree well with the benchmark data while the slight difference in circularity is due to the different filament shapes.Results mesh M3 resolves the principal flow characteristics and captures the predicted shape for a bubble in this regime.The trailing filament behaviour can be seen as a combination of both the 'MooNMD' and 'FreeLIFE' results due to the long thin filament and larger nodule at the edge.It is worth noting that the frequency of reinitialisation has a large effect on the accuracy of the filaments modelled.This is because the cusps have sharp transition zones.In both cases we perform three reinitialisations per time step so a smaller time-step will require a greater number of reinitialisations.

Benchmark: single-phase viscoelastic flow
In order to validate the suitability of our numerical scheme for solving viscoelastic constitutive equations, we will consider the popular benchmark problem of Oldroyd-B flow passed a confined cylinder in a channel.There is a wealth of data available for this case, for example in Fan et al. [46], Alves et al. [47] and Claus et al. [48], all of which are used as a means of comparison.Since this benchmark case is only utilised in order to confirm the accuracy of the viscoelastic constitutive equation (and its associated stabilisation techniques), we consider creeping flow conditions, where: Since this is a single-phase system, all material parameters are no longer phase dependent and the level-set part of the algorithm is not required.We solve the equations in Eq. ( 68) as a mixed system, so we couple our finite elements together into a dual space and approximate velocity u and pressure p simultaneously.The polymeric stress is advanced in time according to Section 3.5.2and we use the SU and DEVSS-G stabilisation techniques detailed in Sections 3.4.2and 3.4.3.The SU parameter is chosen to be γ = 0.05.We enforce inflow conditions on the left, outflow conditions on the right, no-slip conditions on the top wall and cylinder surface and also a symmetry condition on the bottom wall: Fig. 6.Comparison of the predicted centre of mass, rise velocity, circularity and bubble shape for Case 1 in Table 1 with those generated the 'MooNMD' and 'FreeLIFE' pieces of software [42].
A schematic of these boundary conditions and the structure of the meshes used is available in Fig. 8.More detailed information about the simulation can be found in Claus et al. [48].We set β = 0.59, θ = 1 − β and vary the Weissenberg number between 0.1 and 1.We found a stable time-step to be t = 0.005 for all meshes used.The parameter we measure is the drag coefficient on the surface of the cylinder, given below: where n C is the unit normal to the surface of the cylinder.There are some observations made by Claus et al. [48] that we can concur with here.Namely, a very high resolution is required in order to obtain converged values for C D .This is due to the steep polymeric stress gradients adjacent to the boundary of the cylinder, the magnitudes of which only get larger and thus harder to resolve with increasing Weissenberg number.In addition, stabilisation techniques are crucial to obtain converged solutions as there exists a critical Weissenberg number for 0.7 < Wi < 0.8 upon which polymeric forces in the Navier-Stokes equations dominate and our simulations diverge.Claus also notes that simulations diverge with increasing spatial refinement, possibly due to unphysical predictions made by the viscoelastic model or the propagation of numerical errors.We have found the same issue, but like Claus our simulations seem to converge to a final value drag value before diverging.Convergence is evident when looking at Table 4 and Fig 1 obtained with those generated by the 'MooNMD' and 'FreeLIFE' pieces of software [42].

Table 3
The meshes used in Section 4.2 with their associated number of cells, minimum cell length and degrees of freedom (DOF) for each variable.The VMc mesh is not used for approximation, rather just for illustration of the structure of the mesh in Fig. 8 for our range of Weissenberg numbers.Our scheme performs very well for low Weissenberg numbers, but losing some accuracy as a result of increased polymeric activity in the higher Weissenberg number regime.The observed convergence to the values provided in the literature successfully validate the implementation of our viscoelastic flow solver, enabling us to combine this and the results in Section 4.1 together in Section 4.3.

Two-phase viscoelastic flow
In order to fully explore the case of a rising gas bubble in a viscoelastic fluid, we will analyse the effect of changing various material parameters and dimensionless numbers associated with the bulk fluid.These include the polymer concentration, the Weissenberg number, the Eötvös number and the mobility factor in the case of the Giesekus model.We analyse the velocity field and shape of the bubble in order to identify and characterise bubble phenomena.We then look at polymeric stress gradients in order to better determine what causes the phenomena.3.

Table 4
Table comparing different values for the drag on the surface of a cylinder for Weissenberg numbers between 0 and 1.A * next to a drag value indicates the solution did not converge for all time, but converged to an apparent drag value before breaking down.3. We plot our results against the aforementioned studies from the literature.We include a zoomed inset axis to show the convergence properties for the drag values at higher Weissenberg numbers.

Polymer concentration
The polymer concentration c of the viscoelastic fluid contained in the domain 2 is the amount by weight of polymer which has been dissolved in the Newtonian solvent.It is varied in order to determine how the presence of more or less polymer in the fluid affects a bubble's rise, rather than changing the qualitative characteristics of the fluid through parame- ters such as the Weissenberg number Wi.Given a polymer concentration c, the solvent and polymeric viscosities are defined according to the following mass fraction relations: Vahabi et al. [23] employ a polymer concentration of c = 13.286wt.% for Figure 6 in their paper investigating the use of the WC-SPH method on viscoelastic rising bubbles.Vahabi et al. [23] obtained the parameters for their work from an earlier study conducted by Zainali et al. [24], who used another SPH based method, namely the I-SPH method.We explore how bubble shape changes for a standardised c = 13 wt.%(similar to the c = 13.286wt.% chosen in [23] and [24]) and then explore the influence of c on bubble dynamics by performing simulations for c = 6.5 wt.% and c = 26 wt.%.In the dimensionless case, a polymeric concentration of c = 26 wt.% will correspond to a viscosity ratio of β = 0.03698, an incredibly small and almost unphysical value.The effect this has on the momentum equation is a drastic reduction in elliptic Newtonian viscosity generated by The computational domain is the same as that used by Vahabi et al. (Figure 1 [23]), spanned by the co-ordinates [0, 2] × [0, 4] and the centre of the bubble is initially at x c = (1, 1) with radius r = 0.3.Material parameters are given in Table 5.
As we observed in Section 4.1, a Newtonian bubble will rise due to buoyancy forces and deform according to interfacial tension forces.For smaller Reynolds numbers the bubble will deform to become a ellipsoid and for high Reynolds numbers a transition from the ellipsoidal to the spherical cap regime is observed.This is due to inertial forces overcoming the forces due to surface tension on the interface.The same happens when the bulk fluid is viscoelastic, except that it is now the polymeric stress term which heavily influences the shape of the interface as it deforms, also by overcoming the interfacial forces.In Fig. 10 the bubble has only just begun to rise so Newtonian forces are still very present in the fluid, due to polymeric forces taking time to accumulate.The values of the polymeric stress tensor τ p are small (when compared to Figs. 11 and 12) so Newtonian forces are allowed to dominate, causing the bubble to rise similar to Case 1 in Section 4.1.
Fig. 11 captures a crucial moment in the simulation as we observe the viscous and viscoelastic stresses overcome the inherent surface tension forces of the bubble, leading to deformation of the interface.Comparing the τ xx and τ yy components of τ p , between Figs. 10 and 11, the polymeric stress regions become focused at stagnation points on opposite ends of the interface, with values almost twice that in the T = 0.13 case compared to the T = 0.05 case.This displays the transition from Newtonian to viscoelastic rise as the increase in extensional forces due to polymer concentration leads to the formation of a trailing cusp shape.The viscoelastic behaviour continues as time progresses, Fig. 12 shows the profile at T = 0.2 and we can see the trailing cusp has extended even more due to the τ yy stress moving down to the edge of the cusp and increasing in magnitude.
We also observe the onset of the negative wake phenomenon that is ubiquitous for viscoelastic fluids with a high enough polymer concentration c.In Fig. 11, subplot 2A can be seen to be forming the classic conical shaped region of flow in the wake of the bubble, where an inversion in the velocity field directly beneath the bubbles cusp contrasts with an upward flow from buoyancy forces to create dual vortices either side of the wake.This effect is decreased and increased for c = 6.5 wt.% and c = 26 wt.% respectively, where a change in the amount of polymer dissolved in the fluid results in a change of polymeric phenomenological activity.The appearance of a negative wake is not affected by rising polymer concentrations in our case, however its duration and intensity is.We can observe in Fig. 13 that at time T = 0.13, when the negative wake is most active, that the interface is extended in both directions, increasingly more so for rising polymer concentration.As the bubble continues to rise the negative wake begins to disappear until the velocity field is more akin to that of a bubble rising in a Newtonian fluid as the viscoelastic stresses have relaxed.

Parametric study
In the previous section we considered the Oldroyd-B constitutive equation for the purposes of extending on the work of Vahabi et al. [23] and Zainali et al. [24].While the model is suitable for the parameters considered, a major drawback is the presence of a singularity in the extensional viscosity as a function of shear rate.This results in the model giving unphysical predictions for certain values of the shear rate.This can be seen in Fig. 2. The Giesekus model which is used in the succeeding section, has no such singularity and is much more reliable for use in simulations with variable shear rate.We will analyse the change in bubble shape and bubble metrics for varying Eötvös and Weissenberg number, and also the Giesekus mobility factor -a parameter that changes the amount of shear thinning behaviour in the viscoelastic fluid.

Weissenberg number
The Weissenberg number is the dimensionless group defined to be the ratio of the time taken for the fluid to relax back to equilibrium and the characteristic time scale of experiment causing the deformation.It allows us to quantify the relative importance of elastic to viscous forces, so a fluid with Wi 1 behaves much more like a viscous fluid whereas for Wi 1 the response to an applied deformation is much more akin to an elastic solid.In Fig. 14 we examine the effect of varying the Weissenberg number for values Wi = 4, Wi = 8 and Wi = 16 on the final bubble shape and associated bubble metrics.For the bubble shape we immediately notice an increase in Weissenberg number correlates to an increase in the rise velocity and consequently the bubble's final position.The peak in the rise velocity is also delayed for increasing Weissenberg number, due to the viscoelastic stresses taking longer to build up.The same behaviour is observed by Pillapakkam et al. [27] when they increased the dimensional relaxation time parameter λ 1 from 0.1 to 0.2.This same mechanism causes the higher final position, as the bubbles with a larger Weissenberg number spend longer in the semi-Newtonian regime with the dominant forces on the interface being viscous and interfacial, rather than viscoelastic.However, when the viscoelastic stresses have been built up, we see they are of a much larger intensity, indicated by the increased length of the tail of the bubble for larger Weissenberg number.We also observe the thinning of the tail for high Wi, possibly indicating that a reduction in drag on the interface also leads to a higher final position.

Eötvös number
The Eötvös number is the ratio of gravitational forces and forces due to surface tension, with a large Eötvös number resulting in relatively negligible interfacial forces.Fig. 15 displays the effect of Eo on the centre of mass, rise velocity, circularity and shape of the bubble.For small Eötvös numbers, the surface tension forces are strong and the bubble interface is not easily deformed by viscous and viscoelastic stresses.As can be seen in Fig. 15, there is no visible tail for Eo = 5 and it only begins to form properly for Eo = 10.Bubbles in this regime are more akin to those examined in Section 4.1 since they are in the ellipsoidal regime.Venkatesan et al. [28] mention in their paper that a critical Capillary number (Ca = Eo/Re) exists at which the bubble will form a trailing cusp.We observe that the critical Eötvös number is in the range 5 < Eo < 10 in our simulations.We also observe only a minor change in bubble shape between Eo = 35 and Eo = 500.

Mobility factor
Extensional stresses become very large in the wake of a rising bubble.The Giesekus model [49] is able to predict shearthinning and extensional-hardening and therefore is more suitable for predicting the dynamics of bubbles.The mobility factor α controls the shear-thinning behaviour of the bulk fluid.Fig. 16 displays the effect of varying α on the centre of mass, rise velocity, circularity and shape of the bubble.Increasing α increases the degree of shear thinning behaviour in the bulk fluid, which decreases the viscous and viscoelastic stresses on the bubbles interface.As a result, there is less resistance to flow and the bubble is able to rise faster and the final position is higher.Even though the viscoelastic stresses around the tail are smaller for larger α, the tail will still extend because the bubble rises faster, explaining the drop in circularity.

Conclusions and future research
In this paper a stabilised finite element formulation for viscoelastic multiphase flow has been presented in which the conservative level-set method is used to define the interface between the phases.To the best of our knowledge, this is the first implementation of this technique to this class of flows.Material parameters were made to depend on the conservative level-set function to provide a distinction between the phases of flow.A viscoelastic constitutive equation is used to evolve the polymeric stress tensor.A numerical scheme based on the finite-element method was then constructed for the non-dimensional version of the problem.An open-source repository was created which gives access to the code used to generate the results in the paper [36].This provides a useful tool for the multiphase viscoelastic flow community and further exemplifies the efficacy of FEniCS as a finite-element library.
The numerical scheme is applied to several benchmark problems to assess its convergence and stability properties.Benchmark problems with increasing levels of difficulty are considered.First, the conservative level-set method was implemented for the two-phase Newtonian benchmark problem of a 2D rising bubble.Excellent quantitative agreement with the literature for bubble metrics such as centre of mass, rise velocity and circularity was obtained for large density and viscosity ratios.The final shape of the bubble for the inertia dominated case had trailing filaments that were resolved to a very high degree.The second benchmark problem was chosen to assess the viscoelastic component of the scheme and, in particular, its associated stabilisation techniques.The benchmark case of single-phase viscoelastic flow past a cylinder in a channel was considered.Excellent agreement with results from the literature for low Weissenberg numbers was observed.For larger Weissenberg numbers we have good agreement and most importantly, convergence to values from contemporary research.We concur with the conclusion that accurate values for high Weissenberg numbers are heavily dependent on spatial refinement.We also observe the absolute necessity of stabilisation techniques like SU and DEVSS-G for stable simulations in the long time limit.
The final benchmark problem combined the multiphase and viscoelastic components of the scheme to simulate a rising gas bubble in a viscoelastic fluid in 2D.The numerical results were compared with the limited number of results in the literature.The DG and DEVSS-G stabilisation techniques allowed us to model polymer concentrations as large as c = 26wt.%and Weissenberg numbers up to Wi = 16.By examining each component of the polymeric stress tensor we were able to explain the mechanical mechanism that leads to the deformation of the bubble interface in detail.In particular we observed the accumulation of τ yy polymeric stresses in the tail, causing the elongated profile of the bubble.By taking snapshots of the bubble's rise over the course of the simulation we identified when phenomenological activity such as the negative wake formed.Lastly we performed a parametric study where we analysed the various dimensionless quantities associated with Oldroyd-B and Giesekus models such as the Weissenberg and Eötvös numbers.By also modelling the Giesekus mobility factor we gained insight into how the shear-thinning behaviour of the viscoelastic fluid can affect a bubble's rise.
In future work we would like to investigate the conditions under which a jump in the rise velocity for a viscoelastic rising bubble is achieved.This behaviour is unattainable using the planar model implemented in the current framework and requires either a fully 3D model or an axisymmetric formulation to explore this further.In addition, the generalisation to three or more phases and multiphase fluid-structure interaction applications [50] are further priorities.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
was financially supported by EPSRC grant (EP/V040235/1), the Royal Society Newton Advanced Fellowship (NAF/R1/201156) and International Exchanges Award (IEC/NSFC/211143, IES/R2/202095).This research was undertaken using the supercomputing facilities at Cardiff University operated by Advanced Research Computing at Cardiff (ARCCA) on behalf of the Cardiff Supercomputing Facility and the HPC Wales and Supercomputing Wales (SCW) projects.We acknowledge the support of the latter, which is part-funded by the European Regional Development Fund (ERDF) via the Welsh Government.

Fig. 1 .
Fig. 1.Illustration of the conservative level-set function for a two-phase flow set-up.The φ(x) = 0.5 contour is the black line between the two phases.

Fig. 5 .
Fig. 5. Schematic of the rising gas bubble simulation with the level-set function φ plotted with its accompanying colour bar.We indicate the direction of the normal vector n the location of the Dirichlet D and Neumann N boundary conditions, plus the direction of gravity.The mesh on the left hand side is a coarse example of the structure of the mesh we use for the results in Section 4.1.
. 9 which tabulate and display the terminal drag values

Fig. 7 .
Fig. 7.Comparison of the predicted centre of mass, rise relocity, circularity and bubble shape for Case 2 in Table1obtained with those generated by the 'MooNMD' and 'FreeLIFE' pieces of software[42].

Fig. 8 .
Fig. 8. Structured mesh used for the single-phase viscoelastic flow benchmark case in Section 4.2.The mesh used corresponds to VMc in Table3.

Fig. 9 .
Fig. 9. Our calculated drag values C D for their respective Weissenberg numbers Wi, using the meshes listed in Table3.We plot our results against the aforementioned studies from the literature.We include a zoomed inset axis to show the convergence properties for the drag values at higher Weissenberg numbers.

Fig. 14 .
Fig. 14.Influence of Wi on bubble metrics at time T = 0.3.All other parameters are kept the same as the c = 13wt.%polymer concentration case from Section 4.3.1.Here t refers to dimensionless time.

Fig. 15 .
Fig. 15.Influence of Eo on bubble metrics at time T = 0.13.All other parameters are kept the same as the c = 13wt.%polymer concentration case from Section 4.3.1.Here t refers to dimensionless time.

Fig. 16 .
Fig. 16.Influence of α on bubble metrics at time T = 0.25.All other parameters are kept the same as the c = 13wt.%polymer concentration case from Section 4.3.1.Here t refers to dimensionless time.

Table 1
Table of material parameters for test cases taken from Hysing et al. [42]. .

Table 5
Table of material parameters from Vahabi et al. [23] for the associated figures in Section 4.3.1.The star indicates that we vary the viscosity ratio β by varying the polymer concentration c.

β
[24] •D, thus acting to increase the hyperbolic character of this equation.This large increase in viscoelastic stresses necessitates the use of stabilisation techniques, so we employ the DG and DEVSS-G methods described in Sections 3.4.1 and 3.4.3.While the work of Vahabi et al.'s[23]expanded on Zainali et al.[24]by correctly predicting a trailing cusp shape, large polymer concentrations were not considered, since no stabilisation techniques were considered in either paper.