Computational modeling of multiphase viscoelastic and elastoviscoplastic flows
Daulet Izbassarov Marco E. Rosti M. Niazi Ardekani Mohammad Sarabian Sarah Hormozi Luca Brandt Outi Tammisola




1]Daulet Izbassarov*
1]Marco E. Rosti
1]M. Niazi Ardekani
2]Mohammad Sarabian
2]Sarah Hormozi
1]Luca Brandt
1]Outi Tammisola
IZBASSAROV et al
[1]Linné Flow Centre and SeRC (Swedish e-Science Research Centre), KTH Mechanics, SE 100 44 Stockholm, Sweden
[2]Department of Mechanical Engineering, Ohio University, Athens, OH, 45701-2979, USA
*Izbassarov Daulet, Linné Flow Centre and SeRC (Swedish e-Science Research Centre), KTH Mechanics, SE 100 44 Stockholm, Sweden. 
daulet@mech.kth.se

[Summary]In this paper, a three-dimensional numerical solver is developed for suspensions of rigid and soft particles and droplets in viscoelastic and elastoviscoplastic (EVP) fluids. 
The presented algorithm is designed to allow for the first time three-dimensional simulations of inertial and turbulent EVP fluids with a large number particles and droplets. 
This is achieved by combining fast and highly scalable methods such as an FFT-based pressure solver, with the evolution equation for non-Newtonian (including elastoviscoplastic) stresses. 
In this flexible computational framework, the fluid can be modelled by either Oldroyd-B, neo-Hookean, FENE-P, and Saramito EVP models, and the additional equations for the non-Newtonian stresses 
are fully coupled with the flow. The rigid particles are discretized on a moving Lagrangian grid while the flow equations are solved on a fixed Eulerian grid. The solid particles are represented 
by an Immersed Boundary method (IBM) with a computationally efficient direct forcing method allowing simulations of a large numbers of particles. The immersed boundary force is computed at the 
particle surface and then included in the momentum equations as a body force. The droplets and soft particles on the other hand are simulated in a fully Eulerian framework, the former with a 
level-set method to capture the moving interface and the latter with an indicator function. The solver is first validated for various benchmark single-phase and two-phase elastoviscoplastic flow 
problems through comparison with data from the literature. Finally, we present new results on the dynamics of a buoyancy-driven drop in an elastoviscoplastic fluid.










Computational modeling of multiphase viscoelastic and elastoviscoplastic flows
    [
    

==============================================================================








§ INTRODUCTION
 
Elastoviscoplastic (EVP) fluids can be found in geophysical applications, such as mudslides and the tectonic dynamic of the Earth. EVP fluids are also found
in industrial applications such as mining operations, the conversion of biomass into fuel, and the petroleum industry, to name a few. Biological and smart materials can be elastoviscoplastic, making the EVP fluid flows relevant for problems in physiology,
biolocomotion, tissue engineering, and beyond. In most of these applications, we are dealing with multiphase flows <cit.>. Therefore, there is a compelling need to study multiphase flows of EVP fluids and predict their flow dynamics in various situations, including three-dimensional and inertial flows. 
Elastoviscoplastic materials exhibit simultaneously elastic, viscous and plastic properties. At low strains the material exhibits elastic 
deformation, whereas at sufficiently high strains the material experiences irreversible deformation and starts to flow. 
Even conventional  yield-stress  test fluids (such as Carbopol solutions and liquid foams) are shown to exhibit simultaneously elastic, viscous and yield stress behavior. 
Hence, in order to accurately predict the behavior 
of such materials, it is essential to model them as an EVP fluid, rather than an ideal yield-stress fluid (e.g., using the Bingham or the Herschel-Bulkley models).  

There are different types of models that have been proposed for EVP fluids. For instance, Saramito <cit.> proposed a tensorial constitutive law under the Eulerian framework, 
which is based on the combination of the Bingham viscoplastic <cit.> and the Oldroyd viscoelastic models <cit.> in a way 
which satisfies the second law of thermodynamics. This model predicts a Kelvin-Voigt viscoelastic solid (an ideal Hookean solid) response before yielding, when the von Mises criterion 
is not satisfied. Once the strain energy exceeds a threshold value that is specified by the von Mises criterion, the material yields, and the stress field is given by the non-linear viscoelastic constitutive law. 
This model was later improved by the same author <cit.> to account for the shear-thinning behavior of the shear viscosity, and also for the smoothness of the plasticity criterion. 
Moreover, this model is capable of predicting the first normal stress difference along with the yield stress behavior in simple shear flows as a result of combining viscoelasticity and viscoplasticity. 

The prediction of an ideal Hookean solid of Saramito's models <cit.> for the EVP material before yielding causes the model to always predict a zero phase difference between the strain oscillation and the material shear stress, which in turn contributes to vanishing viscous harmonics. This results in an erroneous prediction of zero loss modulus (G^''), 
which is in disagreement with the large amplitude oscillatory shear (LAOS) experiments for identifying and characterizing the properties of the EVP materials <cit.>. 
It was shown by Dimitriou and co-workers <cit.> that for a Carbopol gel (an EVP material), in the limits of small deformation amplitudes, the loss modulus  (G^'') 
is always non-zero and indeed is smaller than the storage modulus (G^') by an order of magnitude. The Isotropic Kinematic Hardening (IKH) idea was then suggested by Dimitriou and co-workers <cit.> and Dimitriou and McKinley <cit.> 
to tackle this problem and to specify the parameters of the models correctly. Based on this concept, the material yield stress builds up and evolves in time together with the flow field, 
where the steady state yield stress is determined via the back stress modulus (a new material parameter) and the deformation of microstructure (a hidden internal dimensionless evolution variable). 
By this method, the energy is allowed to be dissipated, and thus, at small strain amplitudes, it predicts a non-vanishing loss modulus. Recently, a comprehensive IKH constitutive framework has 
been developed to model the thixotropic behavior presents in some practical  EVP materials such as  waxy crude oils <cit.>.

De Souza Mendes <cit.> proposed another constitutive equation for EVP fluids. The basic idea of this model is to modify the classical version of the Oldroyd-B equation, 
where the constant parameters, i.e. the relaxation time (λ_1), the retardation time (λ_2) and the viscosity (η), are replaced with functions of the deformation rate. 
This model reduces to the classical Oldroyd-B equation in the limit of zero shear rate for the unyielded material. Benito and co-workers <cit.> presented another minimal, 
fully tensorial and rheological constitutive equation for EVP fluids. This model predicts the material behaviour as a viscoelastic solid,  capable of deforming substantially before 
yielding, and predicts a viscoelastic fluid after yielding. Moreover, based on the second law of thermodynamics this model has a positive dissipation. 
Recently, Fraggedakis <cit.> performed a systematic comparison of these recently proposed EVP fluid models. The models were tested in simple viscometric flows and against 
available experimental data.

A significant number of numerical studies have analysed purely viscoelastic and purely viscoplastic fluids, but a very limited number accounted for EVP fluids in which neither elastic nor plastic 
effects are negligible. The main reason is that numerical simulations of EVP fluid flows are not a straightforward task due to the inherent non-linearity of the governing equations. 
Nevertheless, numerical simulations can provide quantitative information which is extremely difficult to access by experiments in EVP fluids (for example, velocity fields and stress fields, separated 
into different contributions), and also detailed understanding of the physics of the interaction between particles and droplets in EVP fluids. 

Numerical simulations have already helped to reveal elastic effects in liquid foams and Carbopol. First, Dollet and co-workers <cit.> performed experimental measurements for the flow of liquid foam around a circular obstacle, where they observed an overshoot of the velocity (so-called negative wake) behind the obstacle. Then, Cheddadi <cit.> simulated the flow of an EVP fluid around a circular obstacle by employing Saramito's EVP model <cit.>. The numerical simulation using the EVP model captured the negative wake.  
A purely viscoplastic flow model (Bingham model) on the other hand always predicted fore-aft symmetry and the lack of a negative wake, in contrast with the aforementioned experimental observations. The numerical simulations could hence prove that the negative wake was an elastic effect.
Recently, the loss of the fore-aft symmetry and the formation of the negative wake around a single particle sedimenting in a Carbopol solution  was captured by  transient numerical calculations 
by Fraggedakis and co-workers <cit.> by adopting the EVP tensorial constitutive law of Saramito <cit.>. This was in a quantitative agreement with the experimental 
observations by Holenberg and co-workers for the flow of Carbopol gel <cit.>. The elastic effects on viscoplastic fluid flows have also been addressed in numerical 
simulations of the EVP fluids through an axisymmetric expansion-contraction geometry <cit.> by using the finite element method. It was observed that elasticity alters the shape and 
the position of the yield surface remarkably, and elasticity needs to be included to reach qualitative agreement with experimental observations for the flow of Carbopol aqueous solutions <cit.>. 
Computations in the same geometry have also been performed by implementing the hybrid finite element-finite volume subcell scheme, and combining a regularization approach with the EVP model of 
Saramito <cit.>. Furthermore, Saramito model has been used to simulate the flow of liquid foam in a Taylor-Couette cell <cit.>. 
By adopting the EVP constitutive equation proposed by de Souza Mendes <cit.>, the flow pattern of EVP fluids in a cavity was investigated numerically, and it was demonstrated 
that the elasticity strongly affects the material yield surfaces <cit.>. Recently, De Vita <cit.> numerically investigated the elastoviscoplastic flow through porous media 
by adopting Saramito's model.

The motivation behind this work is to develop an efficient and scalable tool to deal with suspensions of particles and droplets in EVP fluids. In this work, we model an EVP fluid via the constitutive 
law proposed by Saramito <cit.>, which provided excellent results in previous numerical studies of e.g. Carbopol, used in many experiments. 

Multiphase viscoelastic (EV) fluid flows have been studied much more than EVP flows, and indeed some of the results in literature will be used to validate our numerical implementation. To give a few examples of such studies, we list 2D and 3D direct numerical simulations of the dynamics of a rigid single particle <cit.>, two particles <cit.>, multiple particles 
 <cit.>, as well as 
droplets in viscoelastic two-phase flow systems in which one or both phases could be viscoelastic <cit.>, including the case of 
soft particles modeled as a neo-Hookean solid (i.e., a deformable particle is assumed to be a viscoelastic fluid with an infinite relaxation time) <cit.>. 













In the case of a pure visco-plastic (VP) suspending fluid, there is an abundance of computational studies of single and multiple particles <cit.>. Full 3D suspension flows for visco-plastic fluids are time consuming, and thus limited to a few benchmark calculations and lower mesh resolutions 
<cit.>. 
However, 2D suspension flows are feasible <cit.>. 

The key computational challenge  is to resolve the structure of the unyielded regions, 
where the stress is below the yield stress, and to locate the yield surfaces that separate unyielded from yielded regions. Two basic methods are used: regularization and the Augmented Lagrangian 
(AL) approach <cit.>. Regularization tends to be faster, but may still require significant more resources than a Newtonian flow.  AL  approaches, although slower, properly 
resolve the stress fields. This is relevant for resolving important physical features of  suspensions of particles in visco-plastic fluids, e.g. the fact that buoyant particles can be held rigidly in suspension 
<cit.>, the limited influence of multiple particles on each other <cit.>, and the finite arrest time, see Ref. <cit.> for more details.
To overcome these limitations,  researchers have addressed yield stress suspensions from a 
continuum modeling closure perspective, deriving bulk suspension properties that agree with rheological experiments <cit.>. 








The present manuscript is organized as follows. In the next section, the governing equations and the elastoviscoplastic constitutive models for multiphase elastoviscoplastic flows in complex geometries are 
briefly described. In Section <ref>, the numerical methodology is presented. In Section <ref>, the numerical method is validated for various single-phase and two-phase elastoviscoplastic benchmark 
problems, and employed for buoyancy-driven elastoviscoplastic two-phase systems. In this work, we adopt two different IBM schemes to simulate EVP suspension flows which are  modifications and improvements 
of the original IBM scheme proposed by Peskin <cit.>. They are explained in section <ref> in more details. 
Finally some conclusions are drawn in Section <ref>. 



§ MATHEMATICAL FORMULATION


The dynamics of an incompressible flow of two immiscible fluids is governed by the Navier-Stokes equations, written in the non-dimensional form as:  

 
    ∇· u = 0,

 
    ρ( ∂u/∂ t +  u·∇ u) = -∇ p 
    + ∇·μ_s ( ∇ u + ∇ u^T ) + ∇·τ + ρ g +  f,



where u =  u(  x, t ) is the velocity field, p = p (  x, t ) is the pressure field, τ = τ(  x, t ) is an extra stress tensor (defined below) 
and g is a unit vector aligned with gravity or buoyancy. The term f is a body force that is used to numerically impose the boundary conditions at the solid boundaries (particle-laden flow) and at the fluid-fluid interfaces (bubbly flow), as described in sections <ref> and <ref>. Finally, ρ and μ_s are the density and the solvent viscosity of the fluid. 

In the present study, the viscoelastic and elastoviscoplastic effects in the flow are reproduced by the extra stress tensor τ. All the flow models (i.e. the Neo-Hookean, viscoelastic Oldroyd-B, FENE-P and elastoviscoplastic Saramito model) can be expressed with a generic transport equation as 

    (∂ B/∂ t+ u·∇ B-  B·∇ u-∇ u^T· B)
     =a/λ I - ℱ/λ B.

where λ and μ_p are the relaxation time and polymeric viscosity, respectively. The definition of B, ℱ and a used in Eq. <ref> are specified in Table <ref> 
for the different models considered in the present study. In the Neo-Hookean material, G is the shear elastic modulus; this model is analogous to considering the material as a viscoelastic fluid with an
infinite relaxation time λ→∞. In the Saramito model, τ^d is the deviatoric stress tensor and its magnitude is defined as

    |τ^d|= √(1/2τ_ij^dτ_ij^d).

In the FENE-P model, L is the extensibility parameter defined as the ratio of the length of a fully extended polymer dumbbell to its equilibrium length.
From a numerical point of view, therefore, the challenges associated to the solution of equation (<ref>) are similar, independently of the material model considered.





§ NUMERICAL METHOD


In this section, we outline the flow solver which has been previously developed for particle-laden flows <cit.>, for bubbly flows <cit.> and 
for viscoelastic flows <cit.>. The grid is a staggered uniform Cartesian grid in which the velocity nodes are located at the cell faces, while the 
material properties, the pressure and the extra stresses are all located at the cell centers. The flow equations are solved 
using a projection method. The spatial derivatives are approximated 
using second-order central differences, except for the advection terms in Eqs. (<ref>), (<ref>) and (<ref>) where the 
fifth-order WENO or HOUC schemes are applied.




 §.§ Non-Newtonian fluid flow

In a non-Newtonian flow, the transport equation for the extra stress tensor (Eq. <ref>) presents specific challenges. Advection terms such as u·∇ B need a special consideration due to the lack of diffusion terms in the equations. The most common approach is to 
introduce upwinding for the advection terms. However, that approach adds artificial dissipation that can cause the configuration tensor B to lose its positive definiteness, which eventually results 
in a numerical breakdown <cit.>. Min <cit.> tested different spatial discretizations for a polymeric FENE-P fluid and showed that a third-order compact upwind scheme has a 
better performance. Dubief <cit.> have also favored this scheme among the others. In both of these studies, a local artificial diffusion is added where the tensor B experiences a loss of positive definiteness (det(B_ij)<0). This discretization 
scheme works well, but is computationally expensive, because it requires to solve a set of linear equations for each component of the configuration tensor in each direction to calculate the derivatives. In this study we have substituted the compact upwind with an explicit fifth-order WENO scheme <cit.>, a considerably less expensive method that matches 
the performance of the compact scheme as the test case below illustrates; the method has been recently used successfully by Rosti and Brandt for an elastic material <cit.>.

Next, we demonstrate the performance of our method in simulating a non-Newtonian fluid flow. A two-dimensional vortex pair interacting with a wall is simulated in a FENE-P fluid, similarly to Min <cit.>. In this flow, Re and Wi are defined as Re = ρΓ_0/(μ_s + μ_p) = 1800 
and Wi = λΓ_0/δ^2 = 5, where Γ_0 is the initial circulation of the vortex and δ is the initial distance between the vortex pair center and the wall. 
The initial radius of each vortex is 0.145 and the distance between the two centers is set to two radii. The solvent to total viscosity ratio 
β=μ_s/(μ_s + μ_p) is 0.9 and the FENE-P extensibility parameter L^2 is 400. 
Simulations are performed in a domain of size 2δ× 2πδ, with 64 grid cells per δ. Periodic boundary conditions are employed in the x-direction, and no-slip/no penetration boundary conditions are employed in the y direction. 
A time sequence of the vorticity contours is shown in figure <ref>, where the result for a Newtonian flow is also given as a reference. It can be observed that the secondary vortices 
are significantly attenuated in the polymeric flow.

A local artificial diffusion is added to the polymer equations (<ref>) in two instances: if the tensor B experiences a loss of positive definiteness (det(B_ij)<0), and if the trace of the tensor B 
reaches 95% of its maximum (which is L^2). 
It is worth pointing out that in the case shown here, artificial diffusion was added in only a fraction of 0.1% of the grid points. 
Contours of B_xx, B_xy, B_yy and the trace of tensor B, normalized with L^2 are given in figure <ref> at t = 15δ^2/Γ_0. Adding the artificial diffusion to 
only a small fraction of grid points preserves the sharp spatial gradients of the tensor B, as shown in this figure. The required amount of artificial diffusion needs to be tuned for each individual 
simulation as it changes with the relevant parameters of the polymeric flow; e.g. simulating the same test case here with L^2=100 removes any need for local artificial diffusion.








 §.§ Bubbly flow 


Fluid-fluid interfaces are captured by the interface-correction level-set method <cit.>,
and the surface tension force is described by the continuum surface force (CSF) model. The second-order Adams-Bashforth scheme (AB2) is used for the integration of governing equations of an EVP bubbly flow. Note that the AB2 scheme is used to facilitate the implementation of the fast pressure-correction
method developed by Dong and Shen <cit.> and Dodd and Ferrante <cit.>.




  §.§.§ Level-set method

In two-fluid systems, an interface between the phases can be resolved using a fully Eulerian method.  
The body force f due to surface tension, see Eqs. (<ref>), is expressed as:

    f = σκδ(ϕ)  n,

where δ is a regularized delta function and σ is the surface tension.

In this paper, we have adopted a mass-conserving, interface-correction level-set method to capture an interface by a continuous level-set function. 
The level-set function ϕ( x,t) approximates the signed distance from the interface. Hence, ϕ = 0 denotes the interface, ϕ > 0 denotes fluid 1 and ϕ < 0  fluid 2. 
The interface is convected with the local velocity field, 
    ∂ϕ/∂ t +  u·∇ϕ = 0.

To calculate the body force in Eq. <ref>, the unit normal vector, n, and the local mean curvature, κ, can be simply computed as

    n = ∇ϕ/|∇ϕ|,


    κ = -∇· n.


With time, if simply advected, the level set field will no longer equal a signed distance to the interface. It is essential that the signed distance 
property is preserved near the interface, because of the normal and curvature computation. We therefore redistance the level set field every 10-20 time 
steps by solving the Hamilton-Jacobi (reinitialization) equation:

    ∂ϕ/∂𝒯+S(ϕ_0)(|∇ϕ|-1)=0,

where ϕ_0 is the level set field before redistancing, 𝒯 is pseudo-time and S(ϕ_0) is the mollified sign function <cit.>. 
One can observe that if a steady state is reached, then the zero level set contour (interface location) is unaltered, while the level set field has returned 
to a signed distance function. In practice, this equation is iterated only for a few steps towards steady state. 
The level set advection Eq. (<ref>) and reinitialization Eq. (<ref>) are solved using a three-stage total-variation-diminishing (TVD) 
third-order Runge-Kutta scheme <cit.>.

The density, the solvent and the polymeric viscosities, and the relaxation time vary across the fluid interface and are expressed in a mixture form as
 
    ρ   = ρ_1 H(ϕ) + ρ_2 (1-H(ϕ)),   λ = λ_1 H(ϕ) + λ_2(1-H(ϕ)), 
    μ_s    = μ_s,1 H(ϕ) + μ_s,2(1-H(ϕ)),   μ_p = μ_p,1 H(ϕ) + μ_p,2(1-H(ϕ)),

 
where the subscripts 1 and 2 denote the properties of the bulk and suspended fluids, respectively, and H(ϕ) is the regularized Heaviside function defined such that it is zero inside the bubbles 
and unity outside.

A complete description of the level set methodology can be found in Sussman <cit.> and Ge <cit.>, 
and references therein.




  §.§.§ Time integration: Adams-Bashforth scheme


To advance the solution from time level n to n+1, we proceed as follows. First, we advance the level set function and update the density and viscosity fields accordingly. Second, we advance the extra stress tensor and the velocity field in time with the second-order Adams-Bashforth scheme as

    B^n+1 =  B^n+Δ t (3/2 RT_ab^n-1/2 RT_ab^n-1),


    u^* = u^n +Δ t (3/2 RU_ab^n-1/2 RU_ab^n-1),

where we have defined the right-hand sides of the Eqs. (<ref>) RT_ab^n and of eq. (<ref>) as RU_ab^n, with

    RT_ab^n = [1/λ( a  I -
      ℱ B) - 
      ( u·∇ B- B·∇ u-
      ∇ u^T· B) ]^n,


    RU_ab^n =    -∇· ( u u)^n +  g + 1/ρ^n+1(∇·[μ_s^n+1(∇ u^n
    +(∇ u^n)^T)] +∇·τ^n 
      +σκ^n+1δ(ϕ^n+1) n^n+1).

To enforce a divergence-free velocity field, Eq. (<ref>), we proceed by solving the Poisson equation for the pressure <cit.>, 
    ∇·(1/ρ^n+1∇ p^n+1) = 1/Δ t∇· u^*,

and finally, the velocity at the next time level is corrected as

    u^n+1 =  u^* -Δ t/ρ^n+1∇ p^n+1.


In the droplet-laden flow, the pressure Poisson equation is solved in both phases, with unequal densities. Per default, the left hand side of the Poisson equation has variable coefficients. 
In order to utilise an efficient FFT-based pressure solver with constant coefficients <cit.>, we use the following splitting of the pressure  term <cit.>:

    1/ρ^n+1∇ p^n+1→1/ρ_0∇ p^n+1+(1/ρ^n+1-1/ρ_0)∇(2 p^n-p^n-1),

where ρ_0 is the density of the lower density phase (a constant). With this splitting, and after multiplying by ρ_0, the Poisson equation (Eq. <ref>) becomes:

    ∇·∇ p^n+1 =  ∇·[(1-ρ_0/ρ^n+1)∇(2 p^n-p^n-1) ]  + ρ_0/Δ t∇· u^*.

and the velocity correction (Eq. <ref>) transforms to:

    u^n+1 =  u^* -Δ t[1/ρ_0∇ p^n+1+(1/ρ^k-1/ρ_0) ∇(2 p^n-p^n-1) 
    ].





 §.§ Particle-laden flow 

The governing equations of EVP particle-laden flow are integrated in time with a third order Runge-Kutta (RK3) scheme. The RK3 scheme is third order accurate, low storage, and improves the numerical 
stability of the code, allowing for larger time steps. Both rigid and deformable particles are considered here. The rigid particles are included using the immersed boundary method (IBM) that allows us to solve the Navier-Stokes equations on a Cartesian grid despite the presence of particles or complex wall geometries, and has become a popular tool in recent years.
The IBM consists of an extra force, added to the right-hand side of the momentum equations, see Eqs. (<ref>), to mimic boundary conditions, creating virtual boundaries inside the numerical domain. This extra force acts in the vicinity of a solid surface to impose indirectly the 
no-slip/no-penetration (ns/np) boundary condition.

In the case of deformable particles, we use the method described in Rosti and Brandt <cit.>. The solid is an incompressible viscous hyper-elastic material undergoing only the 
isochoric motion, where the hyper-elastic contribution is modeled as a neo-Hookean material, thus satisfying the incompressible Mooney-Rivlin law. To numerically solve the fluid-structure interaction problem, 
we adopt the so called one-continuum formulation <cit.>, where only one set of equations is solved over the whole domain. Thus, at each point of the domain the fluid and solid 
phases are distinguished by the local solid volume fraction ϕ_s, which is equal to 0 in the fluid, 1 in the solid, and between 0 and 1 in the interface cells. The set of equations can be closed 
in a purely Eulerian manner by introducing a transport equation for the volume fraction ϕ_s (the equation is formally the same used in the level-set method, i.e., equation (<ref>)). 
The instantaneous local value of the elastic force is found by solving equation (<ref>), which represents the upper convected derivative of the left Cauchy-Green deformation tensor. The right hand side of equation (<ref>) is identically 
zero for a hyperelastic material <cit.>.




  §.§.§ Time integration: Runge-Kutta scheme


When using the third order Runge-Kutta scheme, the extra stress tensor and the unprojected field are computed by defining RT_rk3^k and RU_rk3^k as

    RT_rk3^k =    ζ^k[1/λ(a I -
      ℱ B) - ( u·∇ B- B·∇ u- 
    ∇ u^T· B) ]^k-1
    
                        +    ξ^k[1/λ(a I -
      ℱ B) - ( u·∇ B- B·∇ u- 
    ∇ u^T· B) ]^k-2,


    RU_rk3^k =    - ζ^k∇· ( u u)^k-1 - ξ^k∇· ( u u)^k-2 
    - 2 α^k/ρ∇ p^k-1 +  α^k/ρ(∇·[μ_s(∇ u+(∇ u)^T)] + ∇·τ)^k
       +  α^k/ρ(∇·[μ_s(∇ u+(∇ u)^T)] + ∇·τ)^k-1 ,

which are the right-hand sides of the Eqs. (<ref>) and (<ref>). Integrating in time yields 

    B^k =  B^k-1+Δ t  RT_rk3^k,


    u^* = u^k-1 +Δ t  RU_rk3^k.

In the previous equations, Δ t is the overall time step from t^n to t^n+1, the superscript ^* is used for the predicted velocity, while the superscript ^k denotes 
the Runge-Kutta substep, with k=0 and k=3 corresponding to times n and n+1. 

The pressure equation that enforces the solenoidal condition on the velocity field is solved via a Fast Poisson Solver

    ∇·∇ψ ^k = ρ/2 α^kΔ t∇· u^*,

and, finally, the pressure and velocity are corrected according to



    p^k   = p^k-1 + ψ ^k, 
     u^k   =  u^* - 2 α^kΔ t/ρ∇ψ ^k,


where ψ is the projection variable, and α, ζ, and ξ are the integration constants, whose values are

    [   α^1 = 4/15   α^2 = 1/15    α^3 = 1/6;             ;   ζ^1 = 8/15   ζ^2 = 5/12    ζ^3 = 3/4;             ;      ξ^1 = 0 ξ^2 = -17/60 ξ^3 = -5/12. ]
 



  §.§.§ Immersed boundary method

The IBM force f cannot be formulated by means of a universal 
equation and therefore IBMs differ in the way  f is computed. The method applied in this solver was first developed by Peskin <cit.> and numerous modifications and improvements have been suggested since 
(for a review, see Ref. <cit.>). In this study we use two different IBM schemes: the volume penalization IBM <cit.> to generate obstacles and complex geometries, and the discrete forcing method for moving particles <cit.> to fully resolve particle suspensions in elastoviscoplastic flows.




  
*Volume penalization IBM 

Kajishima and Breugem <cit.> proposed the volume penalization IBM, where the IBM force f is calculated from the first prediction velocity u^* that is obtained by integrating Eq. (<ref>) in time without considering the IBM force and the pressure correction. The IBM force f and the second prediction velocity u^** are then calculated as follows:

  

    f_ijk  =    ρ α_ijk(u_s - u^*)_ijk/Δ t , 
    u^**_ijk  =    u^*_ijk + Δ t  f_ijk / ρ  ,


where α_ijk is the solid volume fraction in the grid cell with index (i,j,k), varying between 0 (entirely located in the fluid phase) and 1 (entirely located in the solid area) and u_s is the solid interface velocity within this grid cell. Figure <ref> indicates the solid volume fractions (highlighted area) for grid cells around u(i,j) and v(i-1,j-1). Solid boundary in this figure is shown by red dashed line. For non-moving boundaries, u_s is 0 and Eq. <ref> reduce to:

    u^**_ijk = (1 - α_ijk) u^*_ijk .

The second prediction velocity u^** is then used to update the velocities and the pressure following a classical pressure correction scheme <cit.>.


The volume penalization IBM is computationally very efficient, since the solid volume fractions around the velocity points can be calculated at the beginning of the simulation using an accurate method, or they can even be 
extracted directly from a physical sample by magnetic resonance imaging or X-ray computed tomography <cit.>. 



  
*Discrete forcing method for moving particles 

Following the IBM framework, we impose the no-slip/no penetration condition at the particle surfaces (Figure <ref>) by adding an extra force f on the right hand side of the fluid momentum equations (<ref>).
Uhlmann <cit.> developed a computationally efficient IBM to fully resolve particle-laden flows. Breugem <cit.> introduced improvements to this method, making it second order accurate 
in space by applying a multi-direct forcing scheme <cit.> to better approximate the no-slip/no-penetration (ns/np) boundary condition on the surface of the particles and by introducing a slight 
retraction of the grid points on the surface towards the interior. The numerical stability of this method for particle over fluid density ratio near unity was also improved by 
accounting the inertia of the fluid contained within the particles <cit.>. Ardekani <cit.> extended the original method to simulate suspension of spheroidal particles  
with lubrication and contact models for the short-range particle-particle (particle-wall) interactions. 






In this study, we apply the same scheme to fully resolved simulations of particle suspensions in elastoviscoplastic flows.  We apply the IBM force on the predicted velocities u^*, 
which have been obtained as in the single-phase situation. The second prediction velocity u^** is then obtained after the application of the IBM force, and u^** 
substitutes u^* in the pressure correction scheme given in the previous section. The formulation to calculate the second prediction velocity is given here:

  

    U^*_l = ∑_ijku^*_ijkδ_d ( x_ijk - X^k-1_l ) Δ x Δ y Δ z    , 
       F^k -1/2_l = ρ_f  U( X^k-1_l ) - U^*_l/Δ t  , 
       f^  k -1/2_  ijk = ∑_lF^k-1/2_l δ_d ( x_ijk - X^k-1_l ) Δ V_l   , 
       u^** = u^* + Δ t  f^  k -1/2 / ρ_f   ,


where capital letters indicate the variable at a Lagrangian point with index l. In equation (<ref>a), we interpolate the first prediction velocity u^* from the Eulerian grid 
to the Lagrangian points on the surface of the particle, U^*_l, using the regularized Dirac delta function δ_d of Roma <cit.>. This approximated delta function 
essentially replaces the sharp interface with a thin porous shell of width 3Δ x; it preserves the total force and torque on the particle in the interpolation, provided that the Eulerian 
grid is uniform. The IBM force at each Lagrangian point, F^k-1/2_l, is proportional to the difference between the interpolated predicted velocity and the local velocity of the surface 
of the particle (for rigid particles, U_p + ω_p ×r, calculated as shown in the paragraph below). In equation (<ref>c), the IBM forces obtained at the 
Lagrangian points are interpolated back to the Eulerian grid by the same regularized Dirac delta function. In equation (<ref>d), the IBM forces in the Eulerian grid (f^  k -1/2_  ijk) 
are added to the first prediction velocity to obtain the second prediction velocity u^**.

Given the smooth delta function and resolutions typically used, the Eulerian forces obtained from two neighboring Lagrangian points overlap. The multidirect forcing scheme proposed by 
Luo <cit.> is therefore employed to iteratively determine the IBM forces such that the no-slip boundary conditions, U^**≈U, are collectively imposed 
at the Lagrangian grid points. The new second prediction velocity u^** is then obtained by solving the equations above iteratively (typically 3 iterations is enough) using the new 
u^** as u^* at the beginning of the next iteration with equation <ref>b substituted by:

    F^k -1/2_l =  F^k -1/2_l + ρ_f  U( X^k-1_l ) - U^**_l/Δ t  .


The second prediction velocity u^** is then used to update the velocities and the pressure following the procedure described in the previous section.

Taking into account the inertia of the fictitious fluid phase inside the particle volumes, Breugem <cit.> showed that equations for particle motion can be rewritten as:

    ρ_p V_p dU_p/d t   ≈     - ∑_l=1^N_LF_l Δ V_l + 
    ρ_f  d/d t( ∫_V_pud V  )  + ( ρ_p - ρ_f )V_pg + F_c ,  
     [8pt] 
    d( I_p  ω_p) /d t   ≈    - ∑_l=1^N_L(r_l ×F_l ) Δ V_l 
    + ρ_f  d/d t( ∫_V_p(r×u) d V  ) + T_c    .

where U_p and  ω_p are the particle translational and the angular velocity, ρ_p, V_p and I_p are the mass density, volume and moment-of-inertia tensor 
of a particle, and r is the position vector with respect to the center of the particle. The first terms on the right hand side of these equations are the summation of IBM forces and torques 
that act on each Lagrangian point. The second terms account for the translational and angular acceleration of the fluid trapped inside the particle shell. The force term - ρ_f V_pg 
accounts for the hydrostatic pressure with g the gravitational acceleration, and F_c and T_c are the force and torque resulting from particle-particle (particle-wall) 
collisions (see <cit.> for more details). These equations are integrated in time using the Runge-Kutta scheme, as explained in the previous section.









§ VALIDATION 




 §.§ Single-phase flow











  §.§.§ Poiseuille flow of a viscoelastic fluid

The first test case deals with the start-up Poiseuille flow of an Oldroyd-B and FENE-P fluid (Bi=0) in a planar channel. The geometry is a two-dimensional channel bounded by two parallel walls separated by a distance h=L_y, where y denotes the wall-normal direction, and x is the streamwise direction. The fluid is initially at rest and set into motion by applying a sudden constant pressure gradient in the 
streamwise direction. No-slip boundary conditions are applied at the walls. As our method solves the three-dimensional Navier-Stokes equations, we impose periodic boundary conditions in the streamwise and 
spanwise directions to emulate the two-dimensional geometry. The following dimensionless variables are introduced: 
 
    Y = y/h;  τ^* = τ h/u_0 (μ_s+μ_p);   T = t(μ_s+μ_p)/ρ h^2;  V=u/u_0;
      β =μ_s/μ_s+μ_p;  Re =ρ u_0 h/μ_s+μ_p;   
      Wi =λ u_0/h,

where τ^* is a non-dimensional stress, T is a non-dimensional time and the velocity scale is u_0=-h^2dp/dx/8(μ_s+μ_p). The uniform grid has the grid size Δ y=h/180. The time-evolution of the centerline velocity and the wall shear stress is shown for 
the Oldroyd-B fluid case in Fig. <ref>. The velocity and the stress components show oscillating behaviour with overshoots and undershoots before settling down to their fully developed values. 
The steady state profiles for velocity and stress components are shown in Figs. <ref> and <ref> for the Oldroyd-B and FENE-P fluids, respectively. As can be seen from these figures, there is
an excellent agreement between our numerical and existing analytical results.





  §.§.§ Temporally evolving mixing layer of a viscoelastic fluid


The FENE-P model implementation has been validated by simulating a viscoelastic temporally evolving mixing layer flow and by comparing our results with those provided by Min <cit.>. 
We consider the initial velocity field u = 0.5 (tanh y), and trigger the roll-up of the shear layer with a small 2D perturbation. 
The characteristic velocity and length scales are 
Δ u= u_max - u_min and δ = Δ u/(du/dy)_max, respectively. The Reynolds number is fixed at Re=ρδΔ u/μ_s=50 and the Weissenberg number at
Wi=λΔ u /δ=25; moreover, the extensibility L^2 is set to 100, and the solvent viscosity ratio tp β=0.9. The dimensionless time is defined as T=t Δ u/δ. 
The 2D numerical domain has the size 30 δ× 100 δ, discretised by 128 × 384 grid points. Note that the flow configuration and domain are the same used by Min <cit.>. 
Figures <ref>(a - c) show the instantaneous vorticity contours for the Newtonian flow, where we can observe that the initial perturbation grows in time and generates two vortices 
(panel a - T ≈ 20), which subsequently roll up (panel b - T ≈ 60) and eventually merge into one large vortex (panel c - T ≈ 100); the polymeric flow shows a similar behavior. 
The quantitative validation is shown in the bottom panel, where we plot the time history of B_xx in the center of the domain: the symbols represent the literature results, whereas the red line indicates  our numerical data. We find a good agreement of the conformation tensor component time history over the whole vortex merging process.





  §.§.§ Shear flow of an elastoviscoplastic fluid

Next, the method is validated for elastoviscoplastic (EVP) single-phase flows. For this purpose, two test cases are considered. The first case is a simple shear flow.
 Initially, the fluid is at rest and set into 
motion by a constant shear rate  γ̇_0. This test case has a constant dimensionless velocity gradient ∇ u=[[ 0 1; 0 0 ]], 
the Weissenberg number Wi=λγ̇_0=1, the Bingham number Bi=τ_0/μ_0 γ̇_0=1, and the viscosity ratio β=1/9. The time evolution of the stresses is shown in Fig. <ref>(a). The stress components increase as long as the yield criterion is not satisfied; once the criterion is fulfilled (T ≈ 1), the energy starts to dissipate as a result of viscous effects, which is clearly seen in the figure as the slope of the time evolution of the stresses decreases significantly.

The second test case considers the periodic shear flow of an EVP fluid. An oscillatory flow is applied by imposing the shear strain γ_0 sin(ω t), where γ_0 is the strain amplitude and ω is the angular frequency of the oscillation. 
The Weissenberg number is defined as Wi=λω and the Bingham number as Bi= τ_y/(ργ_0 ω) in this case. Computations are performed for two  different Bingham numbers, i.e., Bi=0 and 300. Note that these two values are extreme cases for which the material behaves like a viscoelastic fluid (Bi=0) and like an elastic solid (Bi=300).  The viscoelastic case can be reached at large strain amplitudes γ_0 →∞, whereas the elastic solid behavior is obtained when the amplitude is small γ_0 → 0. 
The other dimensionless parameters of the problem are kept constant at Wi=0.1 and β=0. The evolution of the shear stress component τ_xy is displayed in Fig. <ref>(b) for Bi=0 and 300. As can be seen in these figures, there is a good agreement between our simulation and the analytical solutions, thus indicating an accurate solution of the EVP model equations.




 §.§ Multiphase flow in complex fluids








  §.§.§ Sedimentation of a spherical particle in an elastoviscoplastic fluid










After validating the numerical method for simple viscometric flows, the method is now applied to study the sedimentation of a spherical particle in a channel filled with an EVP 
fluid. This problem exhibits different viscometric flows simultaneously, i.e., biaxial stretching upstream of the particle, shear flow on the sides, and uniaxial extensional flow downstream of it. 
The Saramito model is employed here to facilitate comparison of the present results 
with the numerical results by Fraggedakis <cit.> and with the experimental data by Holenberg <cit.>. 
A single spherical particle of radius R is centered in a domain of size (L_x × L_y × L_z)=(12R × 20R × 12R); a grid of 144×240×144 points is used to discretize the computational 
domain. Periodic boundary conditions are imposed in the x (spanwise) and y (gravity) directions whereas a free slip/no penetration condition is enforced in the z direction.  
The particle starts moving due to the gravity in an otherwise quiescent ambient EVP fluid. Following Fraggedakis <cit.>, the non-dimensional parameters are defined as follows: 
 
    Ar = Δρ^2 g R^3/μ_s+μ_p;    Wi=λΔρ g R/μ_s+μ_p;
       Bn = τ_y/Δρ g R;   ρ^∘=ρ_s/ρ_f,

where Ar, Wi, Bn, and ρ^∘ represent the Archimedes number, the Weissenberg number, the Bingham number, and the density ratio, respectively. In Eq. (<ref>), the density difference is defined 
as Δρ = ρ_f (ρ^∘-1) = ρ_s-ρ_f, where ρ_s and ρ_f are the solid and the fluid densities, respectively. 
The present results are compared with the computational simulations by Fraggedakis <cit.> and the experimental results by Holenberg <cit.>. 
The simulation is performed for Ar=0.03, Wi=1.04, Bn=0.089, ρ^∘=1.38 and β=0.01. Note that the “rough heavy sphere" case is considered in the present study, referring 
to the no-slip boundary condition on the particle. First, a quantitative comparison is conducted based on the steady state settling speed of the particle in the EVP fluid. The terminal velocity 
predicted by the present study, the one in the numerical simulation in Ref. <cit.> and the one observed experimentally <cit.> are reported in Table <ref>. 
As can be seen from the table, the present results could capture the terminal velocity accurately, and indeed the result of the present study deviates from that in other published studies by  less than 4%. 
Next, we present a qualitative comparison of the velocity fields around the spherical particle. The steady state velocity contours normalized with the particle terminal velocity are displayed in 
Fig. <ref>. 
Direct comparisons with the experimental data of <cit.> 
are also included for the yield surface represented by the white markers around the sphere, where the circular and square marks indicate two different experimental series. 
Holenberg <cit.> 
determined the yielded region by means of PIV and PTV techniques, i.e., they defined the yielded region where the velocity magnitude exceeded 10% of the settling velocity. 
To facilitate direct comparison with experimental data, the velocity contours shown in Fig. <ref> are constructed as follows: the distance between the consecutive contour lines is the same and 
equals to 10% of the terminal velocity, starting from 10% to 90% of the velocity. 

Generally, we are in good agreement with the experimental marks <cit.> and 
simulation results shown in Fig. 9 in the work of Fraggedakis et al. <cit.>. Similarly to <cit.>, the current methodology could capture the expected loss of the fore-aft symmetry. On the other hand, a slight discrepancy 
between the present results and those by the aforementioned work <cit.> can be attributed to a different computational box and a lower resolution. Indeed, in the present study a full three-dimensional 
flow is employed, whereas Fraggedakis <cit.> consider an axisymmetric configuration. Moreover, local grid refinement is used in Ref. <cit.>. 
Their very fine grid in the vicinity of the sphere may result in an improved resolution of the yielded region. 
Another reason could be the employment of different boundary conditions in the far-field boundary, the open-boundary condition is used by Fraggedakis <cit.> 
while a periodic boundary condition is employed in the present work.














  §.§.§ Deformable dilute suspension in a shear flow
 



  
*Steady deformation of a neo-Hookean elastic particle in a shear flow 

In this test case, we simulate the flow in a plane Couette geometry. We use a Cartesian uniform mesh in a rectangular box of size 16R × 10R × 16R, 
with 16 grid points per particle radius R. Periodic boundary conditions are imposed in the streamwise (x) and spanwise (z) directions, and the no-slip condition at the walls  (y = -h and y = h), which move in two opposite directions with a constant streamwise velocity ± U = h γ̇. The Reynolds number Re=ργ̇R^2/μ is fixed at 0.1 and the Capillary number Ca=μγ̇/G varied one order of magnitude between 0.05 and 0.5. After the transients die out, the sphere deforms to approximately an ellipsoid, and we therefore characterize these shapes 
by the Taylor parameter (D) and the angle θ. The Taylor deformation parameter is defined as D=(L-B)/(L+B), where L and B are the major and 
minor axis of the equivalent ellipsoid in the middle plane, and θ is the inclination angle with the respect to the streamwise direction.
The steady state values of D and θ are reported in Fig. <ref> for different Ca, and compared with those by 
Villone <cit.>. Similarly to the case of a viscoelastic droplet in a Newtonian medium, deformation as well as the tendency to align with the flow  increases with Ca, 
i.e. with the deformability.  A very good agreement is found between our numerical results and those in the literature. Further validation and details of our implementation can be found in 
<cit.>.





  
*Three-dimensional viscoelastic droplet 

Verhulst <cit.> and Cardinaels <cit.> considered fully three-dimensional shear-driven droplets in which either the droplet or the surrounding fluid is viscoelastic. The Oldroyd-B model is employed in the present study to facilitate a direct comparison with the results of Verhulst <cit.> and Cardinaels <cit.>. 
The spherical droplet of radius R is at the center of the computational domain. Opposite velocities,  V and -V, are enforced on the two walls located at z = 0 and z = H
 to obtain the shear rate γ̇=2V/H. Periodic boundary conditions are imposed in the x (spanwise) and y (streamwise) directions and no-slip conditions at the two walls. 
Following Verhulst <cit.> and Cardinaels <cit.>, the non-dimensional parameters are defined as follows: the Reynolds number Re=ρ_1 γ̇R^2/μ_1, 
the capillary number Ca=R γ̇μ_1/σ, the Weissenberg number Wi=λγ̇, the viscosity ratio k_μ=μ_2/μ_1, the density ratio k_ρ=ρ_2/ρ_1, and the confinement ratio 
χ=2R/H. Alternatively, two Deborah numbers can be defined as De_1=(1-β)Wi/Ca and De_2=(1-β)Wi/(k_μ Ca). 

The results are presented in terms of the Deborah number and dimensionless 
capilary time t γ̇/Ca. The droplet deformation in the y-z plane is measured by the Taylor deformation parameter introduced above.


Following Ramanujan and Pozrikidis <cit.>, the inertia tensor of the drop is used to find the equivalent 
ellipsoid. 

First, we consider the startup dynamics of an Oldroyd-B droplet in a Newtonian medium (VN). 
The viscoelastic spherical droplet is centred in a computational domain of size L_x=H, L_y=2H, L_z=H, 
which is discretised with a resolution of Δ x= Δ y = Δ z = H/192. The simulations are performed at Re=0.05, De_2=1.54, β=0.68, k_ρ=1, k_μ=1.5 and χ=0.25. 
The time evolutions of the Taylor parameter and the angle of inclination for a viscoelastic droplet in a Newtonian fluid at Ca=0.14 and 0.32 are depicted in Fig. <ref>(a) together with the 
numerical results by Verhulst <cit.>. As expected, the drop deformation and alignment with the flow increase with Ca. Also, the time evolution of both the Taylor parameter 
and the inclination angle are in good agreement with the results reported in Ref. <cit.>. 

Next, the dynamics of a Newtonian droplet in an Oldroyd-B fluid (NV) is studied. The resolution is fixed at Δ x= Δ y = Δ z = H/64 and the computational domain is L_x=2H, L_y=4H, L_z=H. 
The computations were performed for Re=0.1, De_1=1, Ca=0.2, β=0.68, k_ρ=1 and k_μ=1.5. The time evolution of the drop deformation parameter and its orientation angle are shown in 
Fig. <ref>(b) for two different confinement ratios: χ=0.46 and χ=0.76. As can be seen in the figure, the confinement ratio increases both the drop deformation and the drop orientation angle. 
The comparison between the present results and those by Cardinaels <cit.> shows good agreement. 




  §.§.§ Buoyancy-driven droplet in viscoelastic and elastoviscoplastic media










Finally, the method is validated for buoyancy-driven (rising) droplets. We start from a Newtonian droplet rising in a Newtonian and viscoelastic fluid.

The Oldroyd-B model is used in the present work to facilitate direct comparison with the results by Prieto <cit.>, Zainali <cit.> 
and Vahabi and Sadeghy 
<cit.>. The fully Newtonian case is a classical benchmark, see e.g., Hysing <cit.>. 
The domain is rectangular with the width L_x=1 and the height L_y=2. A spherical droplet with a radius R is initially placed at the centerline of the channel at a distance of L_y/4 from the lower part of the 
channel. The no-slip boundary conditions are applied at the horizontal walls. It should be noted that Prieto <cit.> used the free-slip boundary conditions on the vertical walls, 
whereas Zainali <cit.> and Vahabi and Sadeghy <cit.> imposed no-slip boundary conditions. The non-dimensional parameters pertaining to this problem are defined as: 
the Reynolds number Re=ρ_1U_g L/μ_1, the Eötvös number Eo=ρ_1 U_g^2 L/σ, the Weissenberg number Wi=λ U_g/L, Bingham number Bi=τ_y/ρ g L, 
the viscosity ratio k_μ=μ_2/μ_1, the density ratio k_ρ=ρ_2/ρ_1 and the confinement ratio χ=2R/L_x. The reference length scale is L=2R, the velocity scale is U_g=√(g L), 
where g is the gravitational constant, and the time scale is L/U_g. 


First, we show a comparison of rising droplets to the computational study of Prieto <cit.> in two cases: (NV) denotes a Newtonian droplet in a viscoelastic medium (Wi=1, β=0.5), 
and (N) denotes a Newtonian droplet in a Newtonian medium (Wi=0, β=1). The other parameters are: Re=35, Eo=10, Bi=0, k_ρ=0.1, k_μ=0.1 and χ=0.5. Figure <ref> 
shows the evolution of the terminal velocity and the steady state shape for a fully Newtonian (N) case and for a Newtonian drop in a viscoelastic medium (NV), both of which are in good agreement with the 
literature results. Note that, in the study of Prieto <cit.>, the microscopic Hooke model was used rather than the Oldroyd-B model considered in the present work; despite of that, very similar results are obtained.

Next, we compare our results in the NV case against the results by Zainali <cit.> and Vahabi and Sadeghy <cit.>. Following these authors, the values of the non-dimensional parameters are Re=1.419, Eo=35.28, Wi=8.083, Bi=0, β=0.07, k_ρ=0.1, k_μ=0.1 and χ=0.3. The droplet interface shapes that we obtained at t=0.13s together with the ones by Zainali <cit.> and Vahabi and Sadeghy <cit.> are depicted in Fig. <ref>. It can be seen 
in the figure that the present result is consistent with the one reported by Vahabi and Sadeghy <cit.>; on the contrary, Zainali <cit.> have not observed the cusped 
trailing edge, which is however a common feature for the case of Newtonian droplet in viscoelastic medium at high polymer concentrations 
<cit.>. 

Finally, some sample simulations are presented for a Newtonian droplet moving in an EVP fluid. The physical properties pertinent to the problem are the same as in Zainali <cit.> 
and Vahabi and Sadeghy <cit.> except for a non-zero Bingham number; indeed, the Bingham number is varied between Bi=0 and 0.1. Fig. <ref> shows shapshots at t=0.13s of the 
streamlines inside the computational domain for the fully Newtonian case (N), and for the EVP fluid with Bi=0, 0.01 and 0.1. In the figure, blue line denotes the solid-fluid boundary defined via the isoline ℱ=0.5, where ℱ is defined in Table <ref>.

Note that all the non-Newtonian cases display a negative wake, and therefore they have four closed streamline zones instead of two zones for the Newtonian droplet.
When increasing the Bingham number Bi, the extent of the yielded region decreases and the solid-fluid boundary approaches the droplet. 
At Bi=0.01, the fluid region occupies most of the domain, but there is a solid region above and below the droplet, as well as two small ellipsoidal regions on both sides of the droplet. 
Finally for Bi=0.1, the solid region occupies almost the whole domain, except for two narrow "caps" at the trailing and leading edges of the droplet. 







§ CONCLUSION

An efficient solver has been presented for the three-dimensional direct numerical simulations of viscoelastic and elastoviscoplastic multiphase flows, expected to allow large-scale simulations also in inertial and turbulent regimes. The solver is general and applicable to non-Newtonian fluids with a dispersed phase which is either rigid or deformable (drops, bubbles and elastic particles). The fluid phases can be chosen to be simple EVP fluids following the model of Saramito<cit.>. The method can be later adapted to more complex EVP models. 

To obtain a stable and accurate solution of the transport equations for the stresses (EVP, elastic or viscoelastic), we use  a fifth-order upwinded WENO scheme for the advection term in the stress model equations. This is found to be very robust and considerably less expensive than the third-order compact upwind scheme suggested in the literature. To avoid numerical breakdown at moderate Weissenberg numbers, a local artificial diffusion can be added. We find that a local diffusion is preferable to the global diffusion which can lead to inaccurate solutions by significantly smearing out the gradients.

The interface between the continuous and dispersed phases is tracked using different approaches for different systems. For the case of deformable viscoelastic particles, we adopt an indicator function based on the local volume fraction. For droplets, we utilise a mass-conserving level set method recently developed by this group, including an accurate computation of the surface tension force based on the local curvature, and a highly efficient and scalable FFT-based pressure solver for density-contrasted flows. The overall solution approach proposed here is independent of the specific interface tracking method. The advantage of these methods is that they are fully Eulerian, efficient, accurate and portable from existing available implementations. For rigid particles, on the other hand, the interface is tracked using an immersed boundary method. In this case, the carrier phase is solved on a fixed Eulerian grid, whereas the interface is represented by a Lagrangian grid following the particle. When comparing to the conventional body fitted grid, the IBM is more simple and versatile for moving rigid bodies.

The method is first validated for single-phase elastoviscoplastic flows including the start up flow in planar channel, temporally evolving mixing layer and simple and oscillating shear flows. Then, it is  applied to the sedimentation of a spherical particle in an EVP fluid, a viscoelastic drop under shear flow and a buoyancy-driven viscoelastic droplet. In all the cases mentioned above, the results obtained with our code are found to be in good agreement with previous results found in the literature. Finally, sample results are presented for a Newtonian droplet rising in an elastoviscoplastic fluid. This, and the behaviour of rigid particle suspensions in EVP fluids, will be interesting topics for future investigations.

The present methodology can also handle multi-body issues. For solid particles, we have a soft-sphere collision model <cit.> and lubrication corrections <cit.>
for short-range particle-particle and particle-wall interactions. In particular, when the gap width between two particles (or particles and wall) reduces to zero,
a soft sphere collision model is activated, to calculate the normal and tangential collision force. We will extend the work on collision models to non-Newtonian fluids in the future. In the level-set method, coalescence takes place automatically. However in some cases, this phenomenon needs to be prevented.
In our previous work <cit.>, a hydrodynamic model was derived for the interaction forces induced by depletion of surfactant micelles. As a future study,
this model could be extended to take into account other effects of surfactants, such as diffusion at the interface and in the bulk fluid.






§ ACKNOWLEDGMENTS

We acknowledge financial support by the Swedish Research Council through grants No. VR2013-5789, No. VR 2014-5001 and VR2017-4809. This work was supported by the European Research Council Grant no. ERC-2013-CoG-616186, TRITOS and by the Microflusa project. This effort receives funding from the European Union Horizon 2020 research and 
innovation programme under Grant Agreement No. 664823.  S.H. acknowledges financial 
support by NSF  (Grant No. CBET-1554044-CAREER), NSF-ERC (Grant No. CBET-1641152 Supplementary CAREER) and ACS PRF (Grant No. 55661-DNI9). The authors acknowledge computer time provided by SNIC 
(Swedish National Infrastructure for Computing) and OSC (Ohio Supercomputer Center).



