J M López-Herrera

Stéphane Popinet

A A Castrejón-Pita

An adaptive solver for viscoelastic incompressible two-phase problems applied to the study of the splashing of weakly viscoelastic droplets

teaching and research institutions in France or abroad, or from public or private research centers.

An adaptive solver for viscoelastic incompressible

two-phase problems applied to the study of the splashing of weakly viscoelastic droplets

Introduction

Using numerical solutions for complex rheologies is nowadays a common predictive tool, since the efficiency of the numerical schemes improves continuously and the computational cost decreases. Three main schemes have been typically used in computational fluid dynamics: Finite Differences (FD), Finite Volume (FV) and finite elements (FE). The presence of interfaces poses additional difficulties. Typical approaches to free surface simulations are the Marker and Cell (MAC), the Volume of Fluid (VoF) and the Level Set (LS) methods. The MAC method has been the reference method for numerous works since the pioneering work of [START_REF] Tomé | A numerical technique for solving unsteady non-newtonian free surface flows[END_REF] [START_REF] Tomé | A numerical technique for solving unsteady non-newtonian free surface flows[END_REF]. Their original version of the MAC scheme is implemented within the framework of the FD method with the advection term approximated using the VONOS scheme [START_REF] Varonos | Development and assessment of a variableorder non-oscillatory scheme for convection term discretization[END_REF]. The original implementation, conceived for simulating Oldroyd-B fluids, has been adapted to solve viscoelastic fluids of finite extensity such as FENE-CR fluids [START_REF] Paulo | Numerical solution of the FENE-CR model in complex flows[END_REF], using the log conformation kernel [START_REF] Martins | A numerical study of the Kernel-conformation transformation for transient viscoelastic fluid flows[END_REF] or the square root kernel [START_REF] Palhares Junior | Numerical study of the square-root conformation tensor formulation for confined and free-surface viscoelastic fluid flows[END_REF]. Other numerical methods, such as the Smoothed-Particle-Hydrodynamics (SPH) method, can also be found in the literature on computational rheology [START_REF] Zainali | Numerical investigation of newtonian and non-newtonian multiphase flows using ISPH method[END_REF][START_REF] Xu | Numerical simulation of 3Dunsteady viscoelastic free surface flows by improved smoothed particle hydrodynamics method[END_REF].

The Finite-Element (FE) method applied to viscoelastic flows goes back to the pioneering work of [START_REF] Nickell | The solution of viscous incompressible jet and free-surface flows using finite-element methods[END_REF][START_REF] Viriyayuthakorn | Finite element simulation of viscoelastic flow[END_REF][START_REF] Guénette | A new mixed finite element method for computing viscoelastic flows[END_REF]. Successful implementations of viscoelastic fluids using FE have recently been conducted [START_REF] Kane | A comparison of four implementations of the log-conformation formulation for viscoelastic fluid flows[END_REF][START_REF] Knechtges | Fully-implicit log-conformation formulation of constitutive laws[END_REF] and is the basis of commercial codes such as Polyflow R . The FE implementation of the log conformation schemes carried out by Hulsen et al. [START_REF] Martien | Flow of viscoelastic fluids past a cylinder at high weissenberg number: Stabilized simulations using matrix logarithms[END_REF] closely follows the publication of the original scheme. The FE implementation of the log conformation performed by Hao & Pan [START_REF] Hao | Simulation for high Weissenberg number[END_REF] is particularly relevant for the present work since the time-split scheme proposed in that work is applied here.

Most of the numerical simulations loose convergence and destabilize when the relaxation parameter, or its dimensionless counterpart, the Weissenberg number, is increased above a threshold value. This behaviour, known as the High-Weissenberg number problem (HWNP), has been a severe hindrance for computational rheology. Fortunately the HWNP can be largely avoided using the technique proposed by Fattal & Kupferman [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF]. These authors proposed formulating the equations in terms of the logarithm of the conformation tensor. Interestingly, this log conformation (kernel) formulation guarantees the positive definiteness of the conformation tensor during the entire simulation. The success of this kernel method has been immediate, and is used in practice as a substitute for the classical approach in computational rheology. The log conformation kernel has been implemented within the FD method [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF], the FE method [START_REF] Martien | Flow of viscoelastic fluids past a cylinder at high weissenberg number: Stabilized simulations using matrix logarithms[END_REF][START_REF] Hao | Simulation for high Weissenberg number[END_REF] and the FV method [START_REF] Afonso | The log-conformation tensor approach in the finite-volume method framework[END_REF]. In the same spirit [START_REF] Balci | Symmetric factorization of the conformation tensor in viscoelastic fluid models[END_REF] proposed using the square root of the conformation tensor to preserve the positive definiteness. Although less common than the log conformation kernel, the square root conformation kernel has been used recently to analyze the lid cavity problem [START_REF] Dalal | Numerical study of driven flows of shear thinning viscoelastic fluids in rectangular cavities[END_REF][START_REF] Palhares Junior | Numerical study of the square-root conformation tensor formulation for confined and free-surface viscoelastic fluid flows[END_REF]. Although other conformation kernels are possible [START_REF] Afonso | The kernel-conformation constitutive laws[END_REF][START_REF] Martins | A numerical study of the Kernel-conformation transformation for transient viscoelastic fluid flows[END_REF], these two seem to be the most accurate.

FV is, at present, the method of reference in CFD (included commercial codes). Several reasons support its popularity. Remarkably, the method is intrinsically conservative, and the simulation of two-phase flows is straightforward since it does not require any special treatment. Among the authors contributing to the development of the FV scheme applied to viscoelastic flows, we can outline Alves and co-workers [START_REF] Alves | The flow of viscoelastic fluids past a cylinder: finite-volume high-resolution methods[END_REF][START_REF] Alves | A convergent and universally bounded interpolation scheme for the treatment of advection[END_REF]. These authors proposed a scheme consisting of a modification of the pressure-velocity coupling SIMPLEC algorithm combined with a new flux limiter for the advection term, CUBISTA, better suited for rheological fluids. Recently it has been shown that a streamfunction-logconformation methodology [START_REF] Comminal | Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/logconformation formulation[END_REF][START_REF] Comminal | Vortex behavior of the Oldroyd-B fluid in the 4-1 planar contraction simulated with the streamfunction-log-conformation formulation[END_REF] can provide stable numerical simulations of flows with very high Weissenberg numbers. Figueiredo et al [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF] have shown that the log-conformation formulation can be used together with the Continuum Surface Force method (CSF) to simulate accurately highly viscoelastic, surface tension dependent, two phase flows. Some implementations are constructed, profiting from existing CFD toolboxes, such as OpenFOAM c [START_REF] Favero | Viscoelastic fluid analysis in internal and in free surface flows using the software OpenFOAM[END_REF][START_REF] Habla | Numerical simulation of viscoelastic two-phase flows using openFOAM R[END_REF][START_REF] Habla | Numerical simulation of the viscoelastic flow in a three-dimensional lid-driven cavity using the log-conformation reformulation in OpenFOAM?[END_REF][START_REF] Pimenta | Stabilization of an open-source finite-volume solver for viscoelastic fluid flows[END_REF]. The efforts of some authors who put their codes at the disposal of the scientific community are worth mentioning. This collaborative spirit allows a continuous improvement of the codes, as those done by Pimenta & Alves [START_REF] Pimenta | Stabilization of an open-source finite-volume solver for viscoelastic fluid flows[END_REF].

Among the most common rheological models, we find the Oldroyd-B [START_REF] Oldroyd | On the formulation of rheological equations of state[END_REF], Giesekus [START_REF] Giesekus | Die elastizität von flüssigkeiten[END_REF], FENE-type [START_REF] Bird | Polymer solution rheology based on a finitely extensible bead-spring chain model[END_REF][START_REF] Chilcott | Creeping flow of dilute polymer solutions past cylinders and spheres[END_REF] or Phan-Thien-Tanner (PTT) [START_REF] Thien | A new constitutive equation derived from network theory[END_REF] models. Each of these models can better suit the particular solvent-polymer solution or melt employed in a given problem. For example, either the Oldroyd-B or the FENE-type seems to fit properly the rheological behaviour of aqueous solutions of polyacrylamide (PAA) [START_REF] Purnode | Polymer solution characterization with the FENE-P model[END_REF][START_REF] Varagnolo | Stretching of viscoelastic drops by steady sliding[END_REF]. Both the FENE-P and FENE-CR models correct the more simple Oldroyd-B model by imposing a maximum stretch that cannot be exceeded (FENE stands for Finitely Extensible Nonlinear Elastic), with the difference between them being the statistical closure used for the restoring force; P denotes the Peterlin's closure [START_REF] Bird | Polymer solution rheology based on a finitely extensible bead-spring chain model[END_REF] and CR follows from the closure proposed by Chilcott & Rallison [START_REF] Chilcott | Creeping flow of dilute polymer solutions past cylinders and spheres[END_REF]. However, numerical simulations seldom match quantitatively the experiments in all of the possible regimes. Note, for example, that numerical simulations, using the Oldroyd-B model, have been employed successfully to explain the origin of the "beads on string" structure appearing in the breakup of weakly viscoelastic droplets [START_REF] Bhat | Formation of beads-on-a-string structures during break-up of viscoelastic filaments[END_REF][START_REF] Turkoz | Axisymmetric simulation of viscoelastic filament thinning with the oldroyd-b model[END_REF], but conversely, overestimate largely the damping factor in slightly vibrating pendant droplets [START_REF] Ponce | On the validity of the jeffreys (oldroyd-b) model to describe the oscillations of a viscoelastic pendant drop[END_REF].

We construct the viscoelastic solver using the free toolbox Basilisk developed by S. Popinet [START_REF] Popinet | Basilisk flow solver and PDE library[END_REF]. Among the different solvers available in Basilisk we can find a library which deals with incompressible fluid problems with a second-order in space time-splitting projection method. Extra forces in the momentum equation can be easily included in the solver in a staggered way to avoid parasitic currents, and facilitate the balance of forces in steady equilibrium situations. The advection term in the momentum equation is computed using the Bell-Colella-Glaz (BCG) second order upwind method [START_REF] Bell | A second-order projection method for the incompressible navier-stokes equations[END_REF]. The VoF method is used for two-phase flows with the advection of the interface performed using the conservative scheme of [START_REF] Weymouth | Conservative volume-of-fluid method for free-surface simulations on cartesian-grids[END_REF]. Surface tension forces are added using Brackbill's CSF procedure [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] in a balanced manner [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. Basilisk also offers tools to easily perform on-the-fly adaptation of the grid depending on the particularities of the flow studied. Adaptation has been used for viscoelastic fluids problems together with FE schemes [START_REF] Saramito | On a modified non-singular log-conformation formulation for Johnson-Segalman viscoelastic fluids[END_REF][START_REF] Jaensson | Simulations of the start-up of shear flow of 2D particle suspensions in viscoelastic fluids: Structure formation and rheology[END_REF]. Saramito [START_REF] Saramito | On a modified non-singular log-conformation formulation for Johnson-Segalman viscoelastic fluids[END_REF] uses an anisotropic auto-adaptive mesh library to search efficient unstructured meshes capable to provide accurate stationary solutions to the lid cavity problem. In the method of Jaensson et al [START_REF] Jaensson | Simulations of the start-up of shear flow of 2D particle suspensions in viscoelastic fluids: Structure formation and rheology[END_REF], the grid moves with the fluid. The mesh tends to become highly distorted and, as a consequence, inaccurate. Jaensson et al tackles the distortion by performing a periodical framing and remeshing as the computation proceeds.

On this platform we have implemented the classic viscoelastic approach in which the advancing equation is written in terms of the stress tensor. We have also implemented kernel conformation approaches, either the log conformation kernel of Fattal & Kupferman, or the square root kernel of Balci et al. In all cases a time-split scheme is used with a calculation of the advection term with the BCG upwind scheme. For the log conformation kernel approach we go further with the time splitting by adopting the scheme of [START_REF] Hao | Simulation for high Weissenberg number[END_REF]. The constitutive model of reference in this work is Oldroyd-B, although for the kernel conformation approaches we have implemented also the FENE-P and the FENE-CR constitutive models for illustrative purpose.

With these implementations we intend to (i) put at the disposal of the scientific community, a validated, ready-to-use, open-source solver using either the log conformation or the square root methodologies that can deal with multiphase flows and fluids of complex rheology; (ii) gain insight on the advantages/drawbacks of the log conformation compared to the square root kernel in the case of two-phase flows; (iii) gain insight on the use of the adaptation of the grid in the resolution of viscoelastic two-phase problems and (iv) report the results on the simulation of the spreading of a weakly viscoelastic fluid after its impact on a flat surface that can be either solid or a liquid layer or bath.

The impact of liquid droplets onto solid surfaces is present in many applications. Most of them search for a control of the coating of the solid by the fluid by managing the dynamics of the impacting droplets. Many investigators have dedicated their efforts to this area of study when the fluid is Newtonian. A thorough review of the state-of-art research on this issue can be found in [START_REF] Josserand | Drop Impact on a Solid Surface[END_REF]. The addition of very small amounts of polymers to a solvent fluid enables a new degree of freedom for this control. In particular, it has been shown that very dilute polymeric solutions inhibit the rebound of droplets over hydrophobic surfaces [START_REF] Bertola | Dynamic wetting of dilute polymer solutions: The case of impacting droplets[END_REF]. In that article the impact dynamics of a droplet of water are compared with that doped with 200 ppm of Polyethylene Oxide (PEO). The spreading stage looks very similar for both Newtonian and viscoelastic fluids. The spreading is dominated by inertia, with negligible viscoelastic forces. Therefore, both droplets reach the same maximum width at the same time. However, the recoiling stage is much slower in the case of the doped droplet. Initially the slowdown of the receding contact line was attributed to the viscoelastic bulk phenomena in the vicinity of the contact line, but direct visualization has shown that the curbing is an interfacial phenomena between the substrate and the drop. The contact line slows down because the polymer molecules are stretched perpendicularly to the contact line as the drop edge sweeps the substrate [START_REF] Bertola | Dynamic wetting of dilute polymer solutions: The case of impacting droplets[END_REF]. Recently [START_REF] Izbassarov | Effects of viscoelasticity on drop impact and spreading on a solid surface[END_REF] [START_REF] Izbassarov | Effects of viscoelasticity on drop impact and spreading on a solid surface[END_REF] and [START_REF] Wang | Impact of viscoelastic droplets[END_REF] [START_REF] Wang | Impact of viscoelastic droplets[END_REF] have studied numerically the spreading and receding of impacting viscoelastic droplets.

[50] use a sharp interface scheme (front tracking) and set the contact angle ad hoc at each computational step with the Kistler correlation. The numerical work described in [START_REF] Wang | Impact of viscoelastic droplets[END_REF] is accomplished using the viscoelastic Giesekus model together with the diffuse-interface Cahn-Hilliard model in which the interfaces are considered as thin transition regions where the interfacial forces are smoothly distributed. The objective of this work was to study the dynamic of the contact line (more precisely the dynamic of the contact angle) when viscoelasticity is present. Recently, [START_REF] Vega | Suppressing prompt splash with polymer additives[END_REF] have studied experimentally the dynamics of the splashing of weakly non-Newtonian drops onto a smooth surface. These authors pay special attention to the transition between different splashing modes (prompt splash, no splash, or corona splash) due to the addition of the polymer. They report that visco-elasticity hinders the development of prompt splashing.

The literature on the splashing of viscoelastic fluids on liquid baths is not as vast as in the case of Newtonian fluids. [START_REF] Pregent | The impact and deformation of a viscoelastic drop at the air-liquid interface[END_REF] studied the splashing of viscoelastic droplets onto either Newtonian or viscoelastic baths. This issue is important, for example, for the formation of capsules or gelled beads. This manuscript is organized as follows. In section 2 the governing equations of an isothermal and incompressible viscoelastic fluid are described. The log and the square root kernels are briefly derived. Details on the numerical schemes are given in section 3. Validation tests of the implemented numerical schemes are performed in section Appendix C. In section 5 we focus on the problem of the splash of weakly viscoelastic droplets.

Governing equations

The equations governing the problem is the set formed by the mass conservation equation,

∇ • u = 0, (1)
and the momentum conservation,

ρ(∂ t u + u • ∇u) = -∇p + ∇ • τ + γκnδ s + ρg (2)
which relates inertia changes to, respectively, the gradient of pressures, fluid internal stresses acting against deformation, surface tension forces and gravitational forces. We denote the density, velocity, pressure and surface tension, curvature by ρ, u, p, γ and κ. δ s stands for the Dirac delta being non-zero at the interface and zero elsewhere. The fluid internal stresses are usually split into the solvent part, τ s , and the polymeric (viscoelastic) contribution τ p ,

τ = τ s + τ p , (3)
while the solvent stress part depends on the deformation tensor as expressed for the usual Newtonian fluid,

τ s = 2µ s D = µ s (∇u + ∇u T ) ,
and the polymeric stress, τ p , takes into account memory effects of the polymers. Several constitutive rheological models are available in the literature with polymeric stresses τ p , which are typically functions f S (•) of the conformation tensor A,

τ p = µ p f S (A)
λ where λ is the relaxtation parameter of the fluid and µ p the polymeric viscosity. The conformation tensor A can be regarded as an internal state variable measuring the molecular deformation of the polymer chains [START_REF] Carreau | Conformation tensor rheological models[END_REF]. The conformation tensor A is assumed to be always symmetric and positive definite, obeying the equation

A = - f R (A) λ (4)
where f R (A) is the relaxation function which is different for each particular constitutive model. denotes the operator upper-convected derivative given by

A = ∂ t A + ∇ • (uA) -A • ∇u -∇u T • A (5) Oldroyd B FENE-P FENE-CR linear PTT f R (A) A -I A 1-tr(A)/L 2 -I A-I 1-tr(A)/L 2 (1 + ε tr(A -I))(A -I) f S (A) A -I A 1-tr(A)/L 2 -I A-I 1-tr(A)/L 2
A -I Table 1: Strain and relaxation functions, f S (A) and f R (A), for some constitutive models [START_REF] Comminal | Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/logconformation formulation[END_REF]. tr(A) stands for the trace of the tensor A.

with ∇u| ij = ∂ i u j . T denotes the "transverse" tensor. In table 1 the expressions of the strain and relaxation functions for some constitutive models are gathered.

Classically, in the case of the Oldroyd-B model, it is usual to skip the use of A by combining Eqs (2) and (4). Then the constitutive equations in terms of the viscoelastic stress tensor, τ p writes,

λτ p + τ p = 2µ p D (6)

The kernel conformation transformation

The numerical resolution of viscoelastic problems often fails to converge when the relaxation parameter, λ, is larger than relatively low values. This instability has been termed in the literature the High-Weissenberg number problem (HWNP), and it has been a major obstacle in computational rheology. Fattal & Kupferman [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF] identified that the instability was caused by a defective modelling of the exponential growths of the stresses. When the instability manifests itself the conformation tensor no longer maintains its property of being definite positive. To tackle the HWNP matrix kernel-transformations of the original conformation tensor have been proposed to enforce at every instant the positive-definite character of the tensor. Two main kernels transformations have been proposed: the log-conformation of Fattal & Kupferman [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF] and the square-root-conformation of Balci et al. [START_REF] Balci | Symmetric factorization of the conformation tensor in viscoelastic fluid models[END_REF].

Log conformation

In this kernel, due to Fattal & Kupferman, rather than advancing the conformation tensor, they suggest to advance in time its logarithm, Ψ = log A. Note that, since A is symmetric and positive-definite, and it is always diagonalizable, then,

A = R Λ R T and Ψ = log A = R log Λ R T ( 7
)
where Λ is the diagonal matrix formed with the eigenvalues and R is the tensor formed by arranging the eigenvectors.

The diagonalization can also be used to decompose the velocity gradient as

(∇u) T = Ω + B + NA -1 (8)
where Ω and N are antisymmetric and B is symmetric, traceless and commutes with A. Using the above decomposition the equation for Ψ is,

∂ t Ψ + u • ∇Ψ -2B -(ΩΨ -ΩΨ) = - e -Ψ λ f R (e Ψ ) (9)
with homogeneous Neumann boundary conditions for Ψ by default. In 2D the decomposition ( 8) is straightforward. In the case of zero polymeric stresses τ p = 0, the elements of the decomposition are, Ω = 0 and

B = 1 2 [(∇u) T + (∇u)].
Otherwise, given the diagonalized conformation tensor

A = R Λ 1 0 0 Λ 2 R T (10)
the velocity gradient is written as

m 11 m 12 m 21 m 22 = R T (∇u) T R (11)
and the elements of the decomposition as

B = R m 11 0 0 m 22 R T , Ω = R 0 ω -ω 0 R T and N = R 0 n -n 0 R T with ω = Λ 2 m 12 + Λ 1 m 21 Λ 2 -Λ 1 and n = m 12 + m 21 Λ -1 2 -Λ -1 1 . (12)
Expressions for the 3D case have been derived in [START_REF] Habla | Numerical simulation of the viscoelastic flow in a three-dimensional lid-driven cavity using the log-conformation reformulation in OpenFOAM?[END_REF]. The square root kernel methodology of Balci et al. [START_REF] Balci | Symmetric factorization of the conformation tensor in viscoelastic fluid models[END_REF] as well as details of its numerical time integration are briefly described in Appendix A while for the classic approach, details are given in Appendix B.

Numerical scheme

We have built the numerical scheme using as a basis the open-source code Basilisk [START_REF] Popinet | Basilisk flow solver and PDE library[END_REF]. Basilisk provides both ready-to-use Finite Volume (FV) solvers for fluid dynamics problems (shallow-water, compressible, incompressible, multiphase...), and an ensemble of useful c-language libraries in order that users can tailor, with a moderate effort, their own specific code.

The incompressible Basilisk solver uses a second-order in space time-splitting projection method. The interface is tracked with a color variable, c(x, t), which represents the volume fraction. c is convected with the fluid,

∂ t c + u • ∇c = 0 . ( 13
)
The above volume fraction equation is solved by successively advecting (sweeping) c along each of the spatial directions, x and y (or r in cylindrical coordinates), using a one-dimensional scheme. As depicted in Fig. 1.a, the onedimensional flux along the sweeping direction is computed from the local linear reconstructed equation, m • x = α, and the face velocities. This one-dimensional net flux must be corrected in case that the one-dimensional velocity field were not divergence-free, i.e, u i-1/2,j = u i+1/2,j in Fig. 1.a. We use the dilation correction proposed in [START_REF] Weymouth | Conservative volume-of-fluid method for free-surface simulations on cartesian-grids[END_REF] which has been proven to be simple, robust and completely volume conservative (if the velocity field is divergence-free). The direction of the first of the one-dimensional sweeps is swapped between x and y in each computational step to avoid a preferred direction of advection.

The surface tension stresses are added to the momentum equation with the CSF method [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] in a balanced manner which avoid parasitic currents [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. The curvature of the interface is computed accurately using the height function approach. In this method the curvature is calculated using the height functions in horizontal or vertical direction, x = h x (y) and y = h y (x), with the curvature κ (say in an almost horizontal interface) given by,

κ = h y 1 + h 2 y .
If the interface is almost vertical, κ can be calculated similarly with x = h x (y) instead of y = h y (x). The method allows to obtain second-order accurate estimates of the curvature. The limits of resolution of the method appear when the size of the cell ∆ is such that κ∆ 1. No special treatment is required in this method for interfacial cells (cells in which the interface is located) next of boundaries and walls. A more detailed description of the method as well as a review of the state of the art in the numerical calculation of surface tension stresses is available in [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Popinet | Numerical Models of Surface Tension[END_REF].

The time stepping of the Navier-Stokes equations is as follows

Step 1: The volume fraction is advanced in time using a conservative, nondiffusive geometric VoF,

c n+1/2 -c n-1/2 ∆t = u n • ∇c n (14)
Step 2: Polymeric stresses are advanced to mid-step n + 1/2, τ n+1/2 p . Step 3: Fluid properties are updated,

θ n+1/2 = θ 1 c n+1/2 + θ 2 (1 -c n+1/2 ) ( 15
)
where θ stands for any property of the fluid; i.e, ρ, µ s , µ p and λ with subscripts 1 and 2 representing the bulk property at each phase.

Step 4: An estimation of the velocity, u * , is calculated by solving

u * -u n ∆t + u n+1/2 • ∇u n+1/2 = 1 ρ n+1/2 -∇p n + ∇ • (2µ n+1/2 D * ) + ∇ • τ n+1/2 p + γκ n+1/2 ∇c n+1/2 (16)
where the advection term is calculated using the Bell-Colella-Glaz second order upwind scheme.

Step 5: The velocity field is projected,

∇ • ∆t ρ n+1/2 ∇p n+1 = ∇ • u * (17)
and updated,

u n+1 = u * -∇p n+1 ∆t/ρ n+1/2 (18)
The time step, ∆t, is determined from two constraints; the stable explicit advection, which implies that the Courant-Friedrich-Levy (CFL) number is below 0.5, and the absence of fake capillary waves which obliges it to have ∆t ≤ (ρh 3 )/(πγ).

Time integration of the polymeric stresses using the log conform kernel

Although we use the log-conformation approach of Fattal & Kupferman we still use as a main variable the polymeric stress tensor, τ p , as proposed by Figueiredo et al. [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF]. Also, in the present scheme we apply the time-split procedure of [START_REF] Hao | Simulation for high Weissenberg number[END_REF] in which Eq. ( 9) is decomposed as

∂ t Ψ + u • ∇Ψ = 0 ( 19
)
∂ t Ψ -2B -(ΩΨ -ΨΩ) = 0 (20)
∂ t Ψ = e -Ψ f R (e Ψ ) λ . (21)
Given the polymeric stresses at time n -1/2, τ n-1/2 p , and the velocity field at instant n, u n a generic time step proceeds as follows:

Step 1: The corresponding conformation tensor at instant n -1/2 is calculated from the relationship,

τ n-1/2 p = λ µ p f S (A n-1/2 )
We assume that the stress function, f S (A), and the relaxation function, f S (A), are linear functions

f S,R (A) = η S,R (ν S,R A -I)
For example, for the FENE-P constitutive model the parameters would be η

S = η R = 1 and ν S = ν R = 1/(1 -T r(A n-1/2 )/L 2 ).
Step 2: The conformation tensor is diagonalised, A = R Λ R T , to obtain its eigenvalues and eigenvectors matrix, Λ n-1/2 and R n-1/2 . Step 3: The log of the conformation tensor is calculated,

Ψ n-1/2 = R log(Λ) R T n-1/2
Step 4: The gradient velocity is decomposed accordingly to Eq. ( 8) to obtain B n and Ω n . Note that for the decomposition we use the eigenvalues and eigenvectors values at instant (n -1/2).

Step 5: The log-conformation tensor is advected using the BCG scheme,

Ψ * = Ψ n-1/2 -∆t∇ • (u n Ψ n )
Step 6: Eq. ( 20) can be integrated explicitly,

Ψ * * = Ψ * + ∆t(2B n + Ω n Ψ n-1/2 -Ψ n-1/2 Ω n )
or, implicitly,

Ψ * * = Ψ * + ∆t(2B n + Ω n Ψ * * -Ψ * * Ω n )
Note that an implicit integration could easily be accomplished given that the resulting equations are linear, and the unknowns at a given point are uncoupled from the unknowns at neighboring points. It would consist in solving N times, once per grid point, a linear system of 3 unknowns (in cartesian 2D; Ψ xx , Ψ yy and Ψ xy ). Our numerical tests on this issue suggest that nothing is gained with the implicit integration. Step 7: The constitutive model Eq. ( 21) is written in terms of the conformation tensor,

∂ t A = - f R (A)
λ and later integrated analytically.

(a) Prior to the analytical integration, the log of the conformation tensor is diagonalised,

Ψ * * = R log(Λ) R T * *
to obtain Λ * * , R * * and the conformation tensor,

A * * = R Λ R T * * . (b) Then, A n+1/2 is calculated with A n+1/2 = A * * e -η R ν R ∆t/λ + (1 -e -η R ν R ∆t/λ ) I ν R
Step 8: Finally,

τ n+1/2 p = µ p λ f R (A n+1/2 ) = µ p λ η R (ν R A n+1/2 -I)

Spatial discretization and the adaptation algorithm

The open source code Basilisk discretizes the computational domain using a structured grid of square finite volumes (termed hereafter cells) that can be either uniform or non-uniform. If a non-uniform grid is preferred, the discretization is arranged hierarchically in a quadtree structure [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries[END_REF] (see Fig. 1.b). In this type of structure, the size of a cell, h, is characterized by its level, , at which is located. Hence, the size of the cells at that level h ∝ 2 -. A prototypical cell of level can be parent of 4 children cells (at the level + 1). The root cell is that corresponding to = 0 from which the rest of the cells at a higher level hang down. A leaf cell is a cell without any child. In the example shown in Fig. 1.b, the grid would be formed by 16 leaf cells with four of them of level = 3, one of level = 1 and the rest of level = 2. All the main variables, including the components of the polymeric stress tensor τ p , are defined at the cell centers. However, the stresses on the right side of the momentum equation ( 16) are computed at the cell faces to avoid any spurious current that could result from the imbalance between pressure and elastic stresses.

(b) (c) (a)
This tree-type grid structure allows a fast and efficient do-loop across the grid nodes. Besides, adding a few constraints in the growth of the tree branches, as for example that the maximum jump of level between neighbouring leaf cells is one, the grid can be refined and coarsened dynamically (adapted) as the simulation proceeds at an affordable computational cost. The adaptation is based on a multi-resolution analysis of selected scalar fields. Consider a control scalar field discretized at grid level , f . This scalar field can be coarsened to the lower level by means of a downsampling operation denominated restriction,

f -1 = restriction(f ) . (22)
This coarser field distribution, f -1 , can be upsampled (or prolongated ) to the original level,

g = prolongation(f -1 ) , (23)
and compared to the original distribution to provide an estimation of the error, ξ = ||f -g ||. Given a particular cell i of level in which the error is ξ i , then that cell will be,

• Refined if ξ i > ζ, • Coarsened if ξ i < 2ζ/3,
• Remain unchanged otherwise.

where ζ is the error threshold set. The prolongation procedure is second-order accurate and involves additional upsampling points in cells contiguous to the finer ones (see Fig. 1.c). A more detailed explanation of the adaptation algorithm can be found in [START_REF] Antoon Van Hooft | Towards adaptive grids for atmospheric boundary-layer simulations[END_REF]. Observe that to fill the new refined and coarsened cells with proper values for each variable can be done with inter/extrapolations that could differ from the prolongation and restriction operators used to decide adaptation regions. In our experience, it is better to use as control adaptation variables the velocity components and the volume fraction. The values of τ p in the newly refined cells are computed with a bilinear interpolation while the coarser values are calculated by averaging. As a result of the hyperbolic nature of the equations for A, boundary conditions ought to be only considered at inflows [START_REF] Van | Mathematical and physical requirements for successful computations with viscoelastic fluid models[END_REF] where we impose by default homogeneous Neumann boundary conditions for tensors, Ψ, b and τ p . However, since in our numerical scheme all the viscoelastic stress components are defined at the center of each cell, care must be taken in order to suitably model the presence of walls and symmetries in the momentum equation [START_REF] Kane | A comparison of four implementations of the log-conformation formulation for viscoelastic fluid flows[END_REF]. Note that the viscoelastic force density applies in our scheme at cell faces and requires to set values at ghost cells since the force density is calculated using central differences. The values at the ghost cells follow the expressions derived in section 3 and 4 of [START_REF] Tomé | Die-swell, splashing drop and a numerical technique for solving the Oldroyd B model for axisymmetric free surface flows[END_REF]. In the case of a rigid wall of orientation n , the normal component, τ p,nn , would be zero. Note that τ p,nn = 0 is only valid for certain constitutive models. For the axisymmetric case, the boundary condition, τ p,θθ = 0 must be added at the wall. Also, on the axis of symmetry the conditions ∂ r τ p,θθ = ∂ r τ p,rr = ∂ r τ p,zz = 0, and τ p,rz = 0 . must be imposed.

Tests

We have performed various tests of the numerical scheme presented in this work in order to verify some aspects of its performance such as the time integration or the correct treatment of the interaction of the viscoelastic fluid with walls and interfaces. Those tests unrelated to the specific splashing problem are gathered in Appendix C.

Splashing of a viscoelastic droplet

This test case is intended to validate the code for axisymmetric two-phase flows in the absence of surface tension. Additionally, some insight on adaptation is gained. The study deals with the time evolution of a viscoelastic Oldroyd-B droplet of density ρ, relaxation parameter λ, solvent and polymeric viscosity µ s and µ p , and diameter D launched from a height H at a velocity U o as sketched in Fig. 2. The surrounding atmosphere is assumed to be dynamically negligible, i.e. ρ a → 0 and µ a → 0. The scaling of the equations of motion will be carried out with the liquid density ρ, the droplet diameter D and the fall velocity U o to give a Froude number F r = gD/U 2 o , a dimensionless height h = H/D, a Reynolds number Re = ρDU o /(µ p + µ s ), a Deborah number De = λU o /D, the ratio of solvent to total viscosity β = µ s /(µ s + µ p ), and the ratio of the outer to inner density and viscosity, ρ r = ρ a /ρ and µ r = µ a /µ o , respectively. This test case has been used by diverse authors with very different schemes [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF][START_REF] Xu | Numerical simulation of 3Dunsteady viscoelastic free surface flows by improved smoothed particle hydrodynamics method[END_REF]. As in the previous work of [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF] the dimensionless parameters were fixed to: F r = 2.26, h = 2, De = 1, Re = 5 and β = 0.1. [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF] do not report values for the outer medium; in the present work we set either µ r and ρ r to 10 -3 . In what follows the dimensional variables are denoted by an asterisk.

We have simulated these tests using the log kernel, the square root kernel and the classic methodology. The computational domain in the present simulations is also shown in Fig. 2. It consists of a square of dimensionless size 2.6 × 2.6. We use axisymmetric equations with the left boundary as the axis of symmetry. The mesh in the simulations is adapted depending on the components of the dimensionless velocity u x , u y and volume fraction f . We have set two ensemble of threshold values; ε th f = ε th ux = ε th uy = 10 -3 (adaptation A1) and ε th f = 10 -3 and ε th ux = ε th uy = 10 -2 (adaptation A2). The simulations performed by [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF] were made with uniform meshes ranging from ∆r = ∆r * /D = ∆z = ∆z * /D = 2.5 × 10 -2 up to ∆r = ∆z = 1.25 × 10 -2 . Since in [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF] negligible differences between meshes are shown, for both adaptation strategies, A1 and A2, the cell widths are comprised between ∆r = ∆z = 2.03 × 10 -2 and ∆r = ∆z = 8.12 × 10 -2 . The maximum timestep has been fixed in all simulations to ∆t = U o ∆t * /D = 10 -3 . Fig. 3 shows the dimensionless width of the droplet, w = W/D versus the dimensionless time t = t * U o /D. We compare our results with the different methodologies against those found in [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF] with the adaptation strategy A1. All three methodologies give very consistent results and are in very good agreement with the results of [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF]. 9) (ii) Eq. (A.1) and (iii) Eq. ( 6). The numerical simulation of [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF] is also shown (continuous red line). The results with the adaptation strategy A2 and the log kernel methodology are also shown.

Splash of weakly viscoelastic drops

In this section we investigate the splash of a viscoelastic drop onto flat substrates. The substrate can be either solid or a viscoelastic liquid film/bath. The properties of the viscoelastic fluid used in the simulations correspond to those of mixtures of pure distilled water with small quantities (around 0.01 wt%) of polymeric solutions of polyacrylamide (PAA), as in the experiments of Vega & Castrejon-Pita [START_REF] Vega | Suppressing prompt splash with polymer additives[END_REF]. Table 2 shows the dependence of the viscoelastic properties, µ p and λ, on the solution concentration. The solvent properties are those of distilled water, µ s = 10 -3 Pa s and ρ = 998 Kg/m 3 . The surface tension is unaffected by the polymeric additives, and is therefore equal to σ = 0.072 N/m. As in subsection 4.1 we use as scaling magnitudes the liquid density ρ, the droplet impact velocity U o and the droplet diameter D (we set D = 3.28 mm as in the experiments of [START_REF] Vega | Suppressing prompt splash with polymer additives[END_REF]). Therefore, a particular splashing is characterized by the following dimensionless quantities:

• A global Reynolds number, Re = ρDU o /(µ s + µ p ).

• A Weber number, W e = ρDU 2 o /σ. Sometimes, in the literature, instead of W e the splashing parameter, K = W e √ Re is used.

• A Deborah number De = λU o /D and a ratio of solvent to total viscosity β = µ s /(µ s + µ p ).

• The ambient to solvent properties ratios, µ r = µ a /µ s and ρ r = ρ a /ρ. • Finally, if the substrate is a liquid film of width L * , its relative depth L = L * /D.

Note that in the above list of parameters the Froude number, F r = U o / √ gD, is absent because it is irrelevant in the splashing phenomena (F r 1) despite the fact that gravity plays a crucial role for accelerating the droplet up to the impact velocity U o . Also other parameters that can be relevant, such as the contact angle or the aspect ratio of the droplet before the impact, are not explored.

The numerical simulation is performed using axisymmetric equations in a square domain similar to the one depicted in Fig. 2. Adaptation is performed at each timestep according to the velocity field and the interface position. The simulations have been carried out with different degrees of grid refinement. Most simulations have been carried out with a grid as fine as 5461 cells per diameter in the adapted region, while far away of that area the grid is coarsened to an equivalent of 21 cells per droplet diameter. Occasionally, for the largest falling velocities, the finest grid reached an equivalent of 10922 cells per diameter. In a few selected cases, the simulations have been performed on parallel machines.

Solid substrate

When the substrate is a solid, the simulation can be started shortly after the impact of the droplet. As shown by [START_REF] Philippi | Drop impact on a solid surface: short-time self-similarity[END_REF], the computed dynamics of the spreading of the droplet, using a slightly truncated landed sphere as initial geometry, is entirely similar to the one obtained while releasing the droplet in the air. While an air dimple can be created when releasing the droplet in air, it does not affect the dynamics of the spreading lamella [START_REF] Philippi | Drop impact on a solid surface: short-time self-similarity[END_REF]. We have selected to initiate the simulations with the the center of the sphere located at a dimensionless distance H/D = (1 -5 × 10 -5 ) above the substrate with a uniform downward dimensionless velocity of the viscoelastic fluid equal to u z = -1. The rest of the variables are set to zero.

To explore the influence of viscoelasticity on the overall dynamics, we focus on the splash of a 1000 ppm solution droplet at an impact velocity U o = 4.09 m/s, that corresponds to W e = 760. For this concentration the other parameters take the following values; Re = 576.33, De = 174.51 and β = 0.043. We impose the ratios values, µ r = 0.018 and ρ r = 0.001. For comparison purposes we also simulate the Newtonian case of a pure solvent (0 ppm).

Fig. 5 and6 show details of the droplet splashing for both 0 and 1000 ppm in concentration. In order to investigate the effect of the wall-fluid interaction we have set a different boundary condition for the volume fraction c, at the wall. The effect of wall-fluid interaction is shown in Fig. 5 where we plot the shape of the lamella at instant t = 0.18. The sliding lamella (red and blue lines) is obtained with the default boundary condition of zero normal derivative: ∂ n c = 0. This condition corresponds to a contact angle of π/2. The same condition, ∂ n c = 0, has been used in Fig. 6 The levitating sheet is obtained by imposing the Dirichlet condition c = 1 i.e. a non-wetting boundary (green and cyan lines). Interestingly, the elastic effects are negligible for these very low polymer concentrations as can be observed in Fig. 5 and6. The mechanism of the splashing is unaffected by the viscoelastic character of the fluid, at least for the very small concentrations cases. It is appealing that the experiments performed by Jung et al. [START_REF] Jung | The role of viscoelasticity in drop impact and spreading for inkjet printing of polymer solution on a wettable surface[END_REF] on the splashing of droplets of solutions of polystyrene in diethyl phthalate over highly wettable solid exhibit the same irrelevance of the polymer concentration in the dynamic of the splashing. Note that both the lamella tip radius, r b , and the radial position of the turning point, r a , are not affected by the viscoelastic stresses in the numerical simulations. Furthermore, the calculated position r a fits well with the analytical Wagner solution obtained using potential theory [START_REF] Philippi | Drop impact on a solid surface: short-time self-similarity[END_REF], similarly to the experiments of Vega & Castrejon-Pita [START_REF] Vega | Suppressing prompt splash with polymer additives[END_REF]. This matching suggests that, in the bulk of the fluid, either viscous and viscoelastic stresses are unimportant during the first stages of impingement.

Viscous and viscoelastic effects are confined to the wall boundary layer and along the contact line. Our numerical results upholds that viscoelasticity could alter the contact line equilibrium that, in turn, affects the dynamics of the lamella. This hypothesis is supported by the experimental results of [START_REF] Jung | The role of viscoelasticity in drop impact and spreading for inkjet printing of polymer solution on a wettable surface[END_REF] where the wettability of the fluids (with and without polymers) is so high that the contact angle is not longer a relevant parameter and the addition of polymers only could be translated into extra bulk elastic stresses. Since the spreading of their droplet is unaffected by the presence of polymers, we might conclude that the bulk elastic stresses have a negligible effect. This is not the general case since the substrate will play a relevant (complicated) role in the spreading and receding stages as the dynamic contact angle will vary during the process [START_REF] Wang | Impact of viscoelastic droplets[END_REF][START_REF] Izbassarov | Effects of viscoelasticity on drop impact and spreading on a solid surface[END_REF]. To summarize, we postulate that viscoelasticity could play a subtle role in the droplet splashing of low-polymer-concentration solutions. The polymer concentration would be too low to promote relevant bulk elastic stresses during the splashing (the elastic stresses would be negligible in the bulk) but sufficient to alter the equilibrium of the contact line and thus cause a distinct behaviour when compared to the Newtonian (zero-concentration) splashing.

Liquid substrate

When the substrate is liquid, the droplet is released at a height equal to H/D = 1.05 setting the dimensionless velocity u z = -1 as an initial condition to all the fluid in the droplet. The other variables are initially zero. The thickness of the film layer has been set to L = 0.3, which seems to be enough to simulate splashing in a deep pool, since simulations done with thicker film layers than L = 0.3 do not show any difference in the mechanism and shape of the splashing. We have simulated splashing with W e ranging from 50 up to 760 that correspond to falling velocities of 1.05 m/s up to 4.09 m/s for a droplet diameter of 3.28 mm, respectively. Fig. 7 depicts the interface evolution of the splashing of a 1000 ppm solution when the impacts are at W e = 50 (right) and W e = 760 (left). The initial instant is t = 0.02 and the timestep is ∆t = 0.06 for both cases. For the most violent impact (W e = 760) an ejecta sheet rapidly develops and creates (after t = 0.14) a well-established corolla. The impact at W e = 50 is completely different. The ejecta sheet is absent and it is replaced by an advancing radial front whose width increases with time (see also the movie in the supplementary material). At a certain time (about t 0.20) the front advances quicker that the fluid bulk and as a consequence the dome formed by the falling droplet becomes surrounded by a kind of toroidal hump instead of a corolla.

Fig. 8 shows the first stages of the splashing for W e = 50 for the pure Newtonian case of 0 ppm and the slightly viscoelastic fluid case of 1000 ppm. In the figure we plot the vorticity distribution ω given by ω = ∂u r ∂z -∂u z ∂r .

(

) 24
Fig. 8 also shows the 2-norm of the conformation tensor Ψ,

||Ψ|| 2 = Ψ 2 rr + Ψ 2 zz + Ψ 2 θθ + 2Ψ 2 rz . ( 25
)
||Ψ|| 2 is used to visualize where the viscoelastic stresses are more intense. As can be seen in the Fig. 8, and the supplementary material, as the drop squeezes the film, the junction front between the drop and film advances and thickens rapidly.

In its advance the front flaps, as a consequence of the vortex shedding, creating a Von Kármán-type vortex street, as was already pointed out by Thoraval et al. [START_REF] Thoraval | Von Kármán vortex street within an impacting drop[END_REF] and confirmed experimentally by [START_REF] Castrejón-Pita | Experimental observation of von kármán vortices during drop impact[END_REF]. At the same time the gas entrapped in the dimple, formed between the droplet and film, rapidly retracts to form a bubble. At the first stages (snapshots t = 0.02 and t = 0.04) no apparent difference exists in the vorticity distribution between the 0 and the 1000 ppm mixtures. However, in subsequent stages it can be observed that the vortex pairs are more distant for the case of a viscoelastic drop (column B) compared to the Newtonian one (column A), since the viscoelastic stresses slightly drag out the shedding of vortices. The evolution of these vortical structures is more interesting. In the case of a Newtonian fluid the vortical structures can only decay by viscous diffusion of the momentum. Since splashing characteristic times are short, and the Reynolds number is large (Re = 3432 for the Newtonian fluid of Fig. 8), the vorticity distribution within the bulk of the liquid is practically the same in snapshot t = 0.06 and subsequent ones. In the case of the mixture of 1000 ppm, the picture is altered by the viscoelastic stresses. Generally speaking, the viscoelastic stresses disrupt this vortical structure as time goes by. Between the spots of positivenegative vorticity, which form the paired vortex, a trail of alternated microvortices appears (shown by the black arrow in the fifth snapshot of column B). Note that, in this case, the spots of vorticity rapidly loose their homogeneity decaying in a turbulent-like mixing. Viscoelastic stresses are concentrated on the fluid surface separating the fluid of the drop from the fluid of the pool, since it is there that larger deformation and strain occur during the splashing process. As can be observed in Fig. 9A, in the lamella, a central core sheet of viscoelastic stresses acts against its spreading and development. In some cases, particularly for violent, high W e number splashes, the viscoelastic stresses tend to bend the incipient lamella, making the first stages of the splashing highly chaotic, as can been seen in sequence 9C. Interestingly, the vortices roll up the viscoelastic stresses resulting in a sort of toroidal spring that delays the advance of the lamella (see Fig. 9A,B andD). These structures are particularly intense when generated around ring bubbles. A sequence of the nucleation of a toroidal spring around a bubble is shown by the green arrow in Fig. 9C. In the first snapshot we can see how the flapping lamella entraps a bag of air by hitting the falling droplet. This bag of air, already has a bubble ring, and is rotating, straining the fluid, and rolling up this strained viscoelastic fluid (second snapshot). Finally, a toroidal spring-like structure is the result. Details of this structure are shown in Fig. 9D.

Effects of Re and β in the first stages of the splashing.

In this section the Weber and Deborah numbers, W e and De are fixed to 300 and 175, respectively, since we intend to explore the role played by Re and β in the first stages of the splashing. We focus on the position of point C that corresponds to the highest position reached by the fluid pertaining to the liquid substrate; the fluid part of the droplet is enclosed by a red line (see the insert in top plot of Fig. 10). In Fig. 10 we explore the role of β on the trajectory of point C for two values of the Reynolds numbers Re = 5 (top) and Re = 500 (bottom). Interestingly, as can be observed in Fig. 10 top, viscoelasticity plays a relevant role for viscous fluids. The imprint caused by the droplet onto the liquid substrate is smaller as the the viscoelastic stresses are more significant (lower values of β). In addition, the splash goes higher for the lower values of β. For a larger Reynolds number, Re = 500, (bottom plot) little effect of viscoelastic stresses is observed in the imprint caused by the droplet (note that the curves depart roughly from the the same position, r C ∼ 0.14). The effect of viscoelaticity on the trajectory of point C is subtle. The trajectory is roughly similar and independent of β but the more important the viscoelasticity (the smaller β), the larger the randomness of the trajectory. Note that the trajectory for β = 0.7 and 0.9 is almost identical and laminar. The effect of a liquid substrate on the spread of the droplet is observed in Fig. 11 where we plot the radial position of point C as a function of time. As can be observed in the bottom plot, viscoelasticity plays a significant role on the size of the dimple imprint at the very beginning of the splashing (roughly for dimensionless times t < 0.07) but once the spreading is launched very little influence on the velocity of the spreading is observed. On the contrary, the overall effect of viscosity is more relevant as can be observed in the top figure. The larger the viscosity of the fluids, the slower the spreading.

Conclusions

In this article we have shown how the time-splitting scheme proposed by Hao & Pan [START_REF] Hao | Simulation for high Weissenberg number[END_REF] can be used together with the classical log conformation tensor of Fattal & Kupferman [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF], or the square-root conformation of Balci et al. [START_REF] Balci | Symmetric factorization of the conformation tensor in viscoelastic fluid models[END_REF], to provide stable numerical simulations of two-phase viscoelastic flows. It is also shown that the time-splitting scheme simplifies the extension of the numerical scheme to different constitutive laws with a moderate effort. Many of the numerical results presented here have been obtained using adaptivity, which can be applied straightforwardly to viscoelastic simulations. The solvers, and most of the tests performed in the present study, are freely available on the Basilisk web page [START_REF] Lopez-Herrera | Viscoelastic log-conform solver and tests[END_REF].

The numerical scheme has been used to investigate numerically the splashing of weakly viscoelastic droplets onto solid flat substrates and pools of the same fluid, taking as reference the experimental conditions of the work of Vega & Castrejon-Pita [START_REF] Vega | Suppressing prompt splash with polymer additives[END_REF]. We observe no difference in the splashing process onto hard substrate between pure solvent droplets and slightly viscoelastic droplets because the viscoelastic bulk effects are negligible for the polymer concentration used. Therefore, we hypothesize that the differences observed by Vega & Castrejon-Pita are due to alterations of the contact line equilibrium because of viscoelasticity and that, in turn, affects the advance of the lamella.

In contrast, the splashing of a slightly viscoelastic droplet onto a pool exhibits a phenomena that has not already been observed in Newtonian fluids. We have observed that the viscoelastic stresses alter the vortex shedding, reported by Thoraval et al. [START_REF] Thoraval | Von Kármán vortex street within an impacting drop[END_REF]. Also, as the splashing proceeds a trail of alternated micro-vortices appears. The viscoelastic stresses are responsible for the disruption of these vortices. The shedding vortices strain the fluid with its rotation, and rolls up this strained viscoelastic fluid to form a kind of toroidal spring. These toroidal springs can nucleate around trapped bubble rings similar to those reported by [START_REF] Thoraval | Von Kármán vortex street within an impacting drop[END_REF]. The size of the dimple imprint on the liquid substrate due to the approaching droplet is affected by viscoelasticity provided the Reynolds number is low enough. The higher the viscoelastic stresses are, the smaller the dimple becomes. Moreover, viscoelasticity causes the a disturbance on the free surface, around the contact line between the drop and the pool, of higher height at low Reynolds number and promotes a random evolution of the lamella at higher Reynolds numbers. The radial velocity of spreading depends solely on the Reynolds number (provided the Deborah and Weber number are fixed).

Appendix B. Numerical scheme for the classic approach

In this scheme we solve Eq. ( 6) by time splitting. The step procedure is as follows.

Step 1: The stress components are advected explicitly with the BCG scheme,

τ * p = τ n-1/2 p + ∆t ∇ • (τ n p u n )
Step 2: The upper convective derivative is solved implicitly, with n = (2k -1)π, α n = 1 + β E n 2 /4 and

G(t) = sinh(θ n t/2) + γ n θ n cosh(θ n t/2) with θ n = α 2 n -E n 2 and γ n = 1 - 2 -β 4 E n 2
where E is the elastic number given by E = λµ o /(ρH 2 ) and β is the ratio of the solvent to total viscosity, β = µ s /µ o = µ s /(µ s +µ p ). In the analytical expression (C.1), the time t and the position y are dimensionless magnitudes. They result, after the time is made dimensionless in λ, t = t * /λ, the y-coordinate with H, y = y * /H and the velocity with the average steady velocity,

u(y, t) = u * (y * , t * ) ū * ∞ being ū * ∞ = - ∆p * ∆x * H 2 3µ o .
where the superscript * denotes the dimensional counterpart.

Using the scaling described above, i.e. H, λ, and ū * ∞ for lengths, times and velocities, respectively, the problem is characterized only by the dimensionless magnitudes E and β being the dimensionless drop of pressure given by,

∆p ∆x = λ ρū * ∞ ∆p * ∆x * = -3E .
Therefore, the numerical simulation domain is a square box of dimensionless size 1 × 1. At the top boundary we set a no-slip condition, while for the bottom symmetry conditions apply. For the left and right boundaries periodic boundary conditions are used for all variables except for the pressure, which is set to 3E at the left side and to 0 at the right side.

We have simulated, using the log kernel approach, the case corresponding to E = 1 and β = 1/9 with three uniform grids with a dimensionless cell size of h = 0.0625 (16 × 16 grid), 0.03125 (32 × 32) and 0.015625 (64 × 64) using a constant time step of value ∆t = 0.001. To use a larger time step compromises the convergence. Subplot A of Fig. C.13 illustrates a comparison between the analytical solution given by Eq. (C.1) with the numerical results obtained with the coarsest grid. As can be observed, the agreement is very good and comparable to similar schemes [START_REF] Habla | Numerical simulation of the viscoelastic flow in a three-dimensional lid-driven cavity using the log-conformation reformulation in OpenFOAM?[END_REF], although the time step in the aforementioned work seems to be smaller. Subplot B illustrates the difference between the theory and the numerical simulation as time proceeds, ε(t) = u(0, t)| theo -u(0, t)| sim for the three grids reported and two different timesteps; ∆t = 10 -5 (continuous line) and ∆t = 10 -3 (dash-point line). The refinement of the grid becomes apparent for ∆t = 10 -3 when the stationary solution is reached. For t = 15 the error with the coarsest grid is 2.82 × 10 -3 dropping to 8.98 × 10 -4 , and to 2.78 × 10 -4 , after each doubling of the spatial resolution. As expected, the error drops with the grid size accordingly to a second-order relation.

The dependence on the viscoelastic model can be observed in Fig. C.14. Subplot A illustrates the temporal evolution of the axial velocity on the axis u(0, t) For two values of the parameter L 2 , L 2 = 10 and L 2 = 1000. For each value of the parameter L calculations has been carried out with the FENE-P and the FENE-CR model. The analytical solution for Oldroyd-B, Eq. (C.1), is also shown. Subplot B illustrates the almost stationary velocity profiles for FENE-P with L 2 = 10, 50 and 1000. As expected the stationary profiles of FENE-CR coincide with the Newtonian parabolic profile. In contrast, the same pressure gradient creates in a FENE-P fluid a larger average velocity (or flowrate) [START_REF] Oliveira | An exact solution for tube and slit flow of a FENE-P fluid[END_REF]. For L 2 → ∞ both FENE-CR and FENE-P coincide with Oldroyd-B. However, the plots in Fig. C.14 show that, in practice, a value L 2 = 1000 suffices.

Appendix C.2. 2D viscoelastic Oldroyd-B droplet immersed in a Couette flow

With this test we wish to validate our scheme using the log kernel methodology when an interface, separating a Newtonian fluid from a viscoelastic one, exists in the presence of surface tension. A sketch of the problem is shown in where the superscript * denotes dimensional variables. The other variables are zero initially. Usually equations are made dimensionless with the outer density, ρ 2 , the droplet radius, a, and the shear rate, γ. With this nondimensionalization, the governing parameters of the problem are: the Weber number W e, the outer Reynolds number, Re, the ratio of viscosities and densities, µ r and ρ r , the Deborah number De, and the ratio of solvent to the total viscosity β given by the following expressions,.

W e = ρ 2 a 3 γ2 σ , Re = ρ 2 a 2 γ µ 2 , µ r = µ 1 µ 2 , ρ r = ρ 1 ρ 2 , De = γ λ and β = µ s µ 1 .
Note that the polymer viscosity is, µ p = µ 1 -µ s and λ is the relaxation parameter. Also, the dimensionless time is, t = γ t * . This problem was first investigated by [START_REF] Chinyoka | Twodimensional study of drop deformation under simple shear for Oldroyd-B liquids[END_REF] and used as a test problem by many others, see for example [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF][START_REF] Khismatullin | Development and implementation of VOF-PROST for 3D viscoelastic liquid-liquid simulations[END_REF][START_REF] Li | Volume preserving immersed boundary methods for two-phase fluid flows[END_REF][START_REF] Zainali | Numerical investigation of newtonian and non-newtonian multiphase flows using ISPH method[END_REF]. In Chinyoka et al. [START_REF] Chinyoka | Twodimensional study of drop deformation under simple shear for Oldroyd-B liquids[END_REF] diverse configurations are explored related to the viscoelastic/Newtonian nature of the outer/inner fluid. Since our objective here is to check how our implementation of the log conformation kernel performs in the presence of a fluid interface, we focus on the configuration with an outer Newtonian fluid surrounded by a viscoelastic where R min and R max are, respectively, the minimum and maximum distance between the interface and the droplet center (the origin in our case). This parameter is also known as the Taylor deformation parameter being denoted by D. very good. In particular, the agreement with the position (t ∼ 0.8) and intensity (equal to approx. 0.0178) of the peak of kinetic energy is excellent. However, the square root kernel has a stationary value (∼ 0.011022) below that obtained with the log kernel (∼ 0.011337) for the M1 grid. When the grid is doubled, i.e. grid M2, the result for the square root kernel increases to ∼ 0.011282 and the log kernel to ∼ 0.011429, closer to the value extracted from [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF] (∼ 0.011572). Adaptation allows the grid to be refined where needed. In the lid cavity problem, as can be observed in Fig. C.18, refinement is located close to the moving wall. Since the velocity is almost established at instant t = 3 the grid distribution shown in Fig. C.18 changes little in later instants. It seems that the log kernel gives a slightly more accurate result than the square root kernel. Interestingly, a similar trend can be observed in Fig. 4.b of [START_REF] Palhares Junior | Numerical study of the square-root conformation tensor formulation for confined and free-surface viscoelastic fluid flows[END_REF]. It is worth that for the calculation of the lid cavity problems, Figueiredo et al. [START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF] report timesteps of 10 -4 , about an order of magnitude larger than ours.

Figure 1 :

 1 Figure 1: (a) Advection scheme of the volume fraction. (b) Quadtree structure (c) Location of the upsampling and downsampling points for the adaptation.

Figure 2 :

 2 Figure 2: Sketch of the problem. The green square denotes the computational domain.

Figure 3 :

 3 Figure3: Time evolution of the dimensionless width of the droplet, w. Numerical simulation using adaptation strategy A1: (i) Eq. (9) (ii) Eq. (A.1) and (iii) Eq. (6). The numerical simulation of[START_REF] Figueiredo | A two-phase solver for complex fluids: Studies of the Weissenberg effect[END_REF] is also shown (continuous red line). The results with the adaptation strategy A2 and the log kernel methodology are also shown.

Figure 4 :

 4 Figure 4: Snapshots at instant t =1.7, 2.2, 3.5 and 5 with adaptation A1 and A2. Each snapshot shows the mesh, the interface position and the spatial distribution of the component of the stress tensor, τ p,θθ

Figure 5 : 2 Figure 6 :

 526 Figure5: Detail of the shape and position of the splash sheet and lamella at t = 0.18 for solutions 0 ppm (Newtonian) and 1000 ppm. In the case of the levitated sheet, the boundary condition at the wall for the volume concentration c, is c = 1. For the sliding lamella the normal derivative of c at the wall is zero: ∂nc = 0. The viscoelastic fluid is modelled using the Oldroyd-B formulation.

19 -

 19

Figure 7 :

 7 Figure 7: Interface evolution in the splashing of a droplet of diameter equal to 3.28 mm of a 1000 ppm solution in a deep pool at falling velocities; left: 4.09 m/s (W e = 760, Re = 1.06, De = 174.51 and β = 7.94 × 10 -5 ) and right: 1.05 m/s (W e = 50, Re = 0.27, De = 44.76 and β = 7.94 × 10 -5). The initial instant is t = 0.02 and the timestep is ∆t = 0.06 in both graphs.

Figure 8 :

 8 Figure 8: Snapshots of the splashing of a droplet with W e = 50. From top to bottom at instants t = 0.02 to t = 0.12 in steps of ∆t = 0.02. Column A shows the vorticity distribution, ω, for 0 ppm (Re = 3432, De = 0 and β = 0). Column B shows the vorticity distribution for 1000 ppm (Re = 0.27, De = 44.76 and β = 7.94 × 10 -5 ). Column C shows the distribution of the 2 norm ||Ψ|| 2 (Eq. (25)) for 1000 ppm.

Figure 9 :

 9 Figure 9: Figures A and B correspond to a 1000 ppm fluid falling with W e = 300 (Re = 0.67, De = 109.64 and β = 7.94×10 -5 ). A and B show ω and ||Ψ|| 2 distributions at instant t = 0.2, respectively. Figure C shows the process of entrapment of a ring bubble and the subsequent roll up of viscoelastic stresses for a 1000 ppm liquid with W e = 760 (Re = 1.06, De = 174.51 and β = 7.94 × 10 -5 ). Figure D shows the details of the structure of the roll up of viscoelastic stresses around ring bubbles.

9 Figure 10 : 9 Figure 11 :

 910911 Figure 10: Trajectory of point C for different values of β. Top: Re = 5; Bottom: Re = 500. The insert in the top plot shows the imprint caused in the liquid substrate (blue line) by the droplet (red line).

Figure C. 12 :

 12 Figure C.12: Sketch of the transient onset planar Poiseuille flow for a viscoelastic fluid. The simulation domain is depicted in green.

1 + λ ∆t τ n+1/ 2 p- 2 p+ τ n+1/ 2 pFigure C. 13 :

 122213 Figure C.13: Subplot (A): Time evolution of axial velocity u(0, t). Comparison between the analytical solution given by Eq. (C.1) and the numerical solution with a uniform grid of 16 × 16 and a timestep ∆t = 10 -3 . Subplot (B): Error evolution for three different grids. The continuous line is obtained with a timestep ∆t = 10 -5 while the dash-point line is obtained with ∆t = 10 -3 .

Figure C. 14 :

 14 Figure C.14: Subplot A: Time evolution of the velocity on the axis for FENE-P and FENE-CR. (L 2 = 10 and 1000). The Oldroyd-B analytical solution given by Eq.(C.1) is also shown. Subplot B: Stationary velocity profile for a FENE-P fluid for L 2 = 10, 50 and 1000. The stationary solution is supposed to be reached at t = 10.

Figure C. 15 :

 15 Figure C.15: Sketch of the benchmark problem of a viscoelastic 2D droplet immersed in a Newtonian fluid undergoing a Couette flow.

Figure C. 16 :

 16 Figure C.16: Plot (A): Deformation Φ versus the dimensionless time t. Plot (B): deformation of the droplet at time t = 10. Dimensionless values of the tested case: Re = 0.3, W e = 0.18, µr = ρr = 1, De = 0.4 and β = 0.5. In both plots the results of the proposed scheme are shown with a continuous black line labelled with the name 'Basilisk'. The numerical simulation has been performed with the log kernel methodology. The open circles correspond to the results of [27].

Figure C. 17 :

 17 Figure C.17: Time evolution of the dimensionless kinetic energy for W i = 1, β = 0.5 and Re = 0.01. The simulations have been performed with two uniform grids; M1 and M2 (64×64 and 128×128, respectively) and an adapted grid. Numerical results with (i) the log-conformation kernel and (ii) the square root kernel are shown. Results of Fattal & Kupferman (2005) are also shown (red circles). Insert: Sketch of the lid cavity problem.

Figure C. 18 :

 18 Figure C.18: Spatial distribution of Ψxy at instant t = 3. The adapted grid at that instant is also shown.

Table 2 :

 2 Polymeric viscosity µp, and relaxation time λ, in S.I units for different diluted PAA-water solutions[START_REF] Vega | Suppressing prompt splash with polymer additives[END_REF].

Acknowledgments

This work has been supported by the Spanish Ministry of Economy Grant DPI2013-46485. J. L-H wishes to thank M.-J. Thoraval for his guidance on the analysis of drop splashing. AACP was funded through an EPSRC-UK grant (EP/P024173/1) and a Royal Society University Research Fellowship.

Appendix A.1. Equations Balci et al. [START_REF] Balci | Symmetric factorization of the conformation tensor in viscoelastic fluid models[END_REF] propose to formulate the constitutive differential models in terms of the (unique) positive symmetric square root b(x, t) of the conformation tensor A(x, t), A = bb T that, substituted in Eq. ( 4), results in the following time evolution equation for b,

where a is an antisymmetric tensor in which off-axis values result from the enforcement of the symmetric character of b. In 2D this would be

Balci et al. [START_REF] Balci | Symmetric factorization of the conformation tensor in viscoelastic fluid models[END_REF] provide expressions for the 3D case.

Appendix A.1.1. Numerical scheme As for the case of the log kernel the numerical scheme is a time splitting procedure of Eq. (A.1). Therefore, a time step can be decomposed in the following substeps.

Step 1: The square root tensor is advected explicitly with the BCG scheme,

Step 2: The rest of Eq. (A.1) is linearized and solved implicitly. Assuming a linear relationship for the relaxation function, the system to be solved would be,

Step 3: Finally the polymeric stress is computed from b n+1/2

. Transient planar Poiseuille flow for a viscoelastic fluid

The problem is sketched on Fig. C.12. A viscoelastic fluid of density ρ, solvent and polymeric viscosity, µ s and µ p , and relaxation time λ, is trapped in the gap of width 2H formed by two parallel infinite plates. The fluid, initially at rest, is set in motion by the sudden application of a constant pressure gradient. The steady planar parabolic Poiseuille flow is reached after a transient period if the viscoelastic fluid is an Oldroyd-B or a FENE-CR. For a FENE-P the profile departs slightly from a strict parabola [START_REF] Oliveira | An exact solution for tube and slit flow of a FENE-P fluid[END_REF]. The transient for a Newtonian fluid is characterized by an exponential increase of the axial velocity. However, in the case of a viscoelastic fluid its elastic nature gives a different behaviour, since an oscillation is superposed to the exponential increase. For an Oldroyd-B fluid, an analytical solution due to Waters & King [START_REF] Waters | The unsteady flow of an elastico-viscous liquid in a straight pipe of circular cross section[END_REF] is available, a square cavity of size L, bounded by walls, except on the top side where a time-dependent tangential velocity is imposed. Using as scaling magnitudes the density ρ, the largest stationary velocity U o and the width of the cavity L, we form a Weissenberg number, W i, a Reynolds number, Re, and a solvent viscosity ratio, β, given by

The standard problem relies on the following regularized dimensionless parabolic profile for the top lid

where x = x * /L, t = U o t * /L and u x = u * x /U o are the corresponding dimensionless variables. The remaining cavity walls are stationary and the no-slip boundary condition is imposed on the four walls. We assume the Stokes limit for the momentum equation. In the simulation, we have set Re = 0.01, β = 0.5 and W i = 1. This test case has become a classical benchmark problem in computational rheology since the HPWN manifests itself with these values of the dimensionless parameters. In Table 1 of [START_REF] Sousa | Lid-driven cavity flow of viscoelastic liquids[END_REF] are gathered previous numerical studies concerned with a lid-driven cavity flow of constant viscosity viscoelastic fluids. We have solved this test case with a uniform grid of 64 × 64 (grid M1) and with a grid of 128 × 128 (grid M2) equivalent to a level = 6 and 7, respectively. The maximum timestep for the M1 grid is ∆t = 5 × 10 -5 , while for M2 we had to set ∆t = 10 -5 . Numerical simulations with an adapted grid have also been carried out. The adaptation is applied every 50 timesteps by controlling the error on the components of the dimensionless velocity. The threshold value for both components is 5 × 10 -4 with the maximum and minimum levels of refinement/coarsening, = 7 and 5, respectively.


Claire Marie

Soleil Dufour

Neville Pillay

Nico Avenant

Johan Watson

Etienne Loire

Guila Ganem
email: guila.ganem@umontpellier.fr

Claire M S Dufour
email: clairems.dufour@gmail.com

Ganem

Orcid Dufour

Ganem Orcid

Guila

Habitat characteristics and species interference influence space use and nest-site occupancy: implications for social variation in two sister species

Keywords: habitat selection, home range overlap, interspecific competition, nest-site sharing, 48 radiotracking, Rhabdomys, secondary contact, social groups 49 50

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.