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

We propose an adaptive numerical solver for the study of viscoelastic 2D two-phase flows using the volume-of-fluid method. The scheme uses the robust log conformation tensor technique of Fattal&Kupferman (2004,2005} combined with the time-split scheme proposed by Hao&Pan (2007}. The use of this time-split scheme has been proven to increase the stability of the numerical computation of two-phase flows. We show that the adaptive computational technique can be used to simulate viscoelastic flows efficiently. The solver is coded using the open-source libraries provided by the \basilisk \cite{Basilisk} platform. In particular, the method is implemented for Oldroyd-B type viscoelastic fluids and related models (FENE-P and FENE-CR). The numerical scheme is then used to study the splashing of weakly viscoelastic drops. The solvers and tests of this work are freely available on the Basilisk web site


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. Typically, three main schemes have been 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 Tomé et al (1996) [6]. Their original implementation of the MAC scheme is implemented within the framework of the FD method with the advection term approximated using the VONOS scheme [7]. The original implementation, conceived for simulating Oldroyd-B fluids, has been adapted to solve viscoelastic fluid of finite extensity as FENE-CR fluids [8], using the log conformation kernel [9] or the square root kernel [10]. Other numerical methods, such Giesekus [33], FENE-type [34,35] or Phan-Thien-Tanner (PTT) [36] models. Each of these models can better suit the particular solvent-polymer solution or melt employed in a particular problem. For example, either the Oldroyd-B or the FENE-type seems to fit properly the rheological behaviour of aqueous solutions of polyacrylamide (PAA) [37,38]. 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 [34] and CR follows from the closure proposed by Chilcott & Rallison [35]. 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 [39], but conversely, overestimate, largely, the damping factor in slightly vibrating pendant droplets [40].
We construct the viscoelastic solver using the free toolbox Basilisk developed by S. Popinet [4]. Among the different solvers available in Basilisk we can find a library which deals with incompressible fluid problems with a second order in a 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 [41]. The VoF method is used for two-phase flows with the advection of the interface performed using the conservative scheme of [42]. Surface tension forces are added using the Brackbill's CSF procedure [43] in a balanced manner [44]. Basilisk also offers tools to easily perform an 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 [45,46]. Saramito [45] 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. [46], the grid moves with the fluid. The mesh tends to become highly distorted and, in consequence, inaccurate. Jaensson er al. tackles the distortion by performing periodically a framing and a remeshing as the computation proceed.
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 Hao & Pan (2007). 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 propose.
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 adaption of 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 impacted 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 [47]. 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 [48]. 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 [48]. Recently Izbassarov & Muradoglu (2016) [49] and Wang et al. (2017) [50] has afforded numerically the study of the spreading and receding of impacting viscoelastic droplets. The authors of [49] 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 [50] 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, [51] have studied experimentally the dynamics of the splashing of slightly non-Newtonian drops onto a smooth surface. These authors pay special attention to the change on the map of the splashing modes (prompt splash, no splash or corona splash) due to the added polymer. They report that visco-elasticity hinders the development of prompt splashing.
The literature on the splashing on liquid baths is not as vast, as in the case of Newtonian fluids. [52] 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, and the momentum conservation, which relates inertia changes to, respectively, the gradient of pressures, fluid internal stresses acting against deformation, surface tension forces and, eventually, gravitational forces. We denote the density, velocity, pressure and surface tension, curvature(which is normal to the interface n) by ρ, u, p, γ and κ. δ s stands for the Dirac delta being one at the interface and zero elsewhere. The fluid internal stresses are usually split into the solvent part, τ s , and the polymeric (viscoelastic) contribution τ p , while the solvent stress part depends on the deformation tensor as expressed for a usual Newtonian fluid, And the polymeric stress, τ p , takes into account memory effects of the polymers. Several constitutive rheological models are available in the literature with their polymeric stresses τ p , which are typically functions f S (·) of the conformation tensor 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 [53]. The conformation tensor A is assumed to be always symmetric and positive definite, obeying the equation where f R (A) is the relaxation function which is different for each particular constitutive model. denotes the operator upper-convected derivative given by Table 1: Strain and relaxation functions, f S (A) and f R (A), for some constitutive models [25]. 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,

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 [1,2] 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 [1,2] and the square-root-conformation of Balci et al. [20].

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, 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 where Ω and N are antisymmetric and B is symmetric, traceless and commutes with A. Using the above decomposition the equation for Ψ is, 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 the velocity gradient is written as and the elements of the decomposition as Expressions for the 3D case have been derived in [30]. The square root kernel methodology of Balci et al. [20] as well as details of its numerical time integration are briefly described in Appendix A while for the classic approach, details are described in Appendix B.

Numerical scheme
We have built the numerical scheme using as a basis the open-source code Basilisk [4]. 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, 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 it is depicted in figure  1.a, the one-dimensional 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 figure 1.a. We use the dilation correction proposed in [42] which it has been proved to be simple, robust 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 preferred direction of advection.
The surface tension stresses are added to the momentum equation with the CSF method [43] in a balanced manner which avoid parasitic currents [44]. 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), being the curvature κ (say in an almost horizontal interface) given by, 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 revision of the state of the art in the numerical calculation of surface tension stresses is available in [44,54].
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, Step 2: Polymeric stresses are advanced to mid-step n + 1/2, τ n+1/2 p .
Step 3: Fluid properties are updated, 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 where the advection term is calculated using the Bell-Colella-Glaz second order upwind scheme.
Step 5: The velocity field is projected, and updated, 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 this has been proposed by Figueiredo et al. [27]. Also, in the present scheme we apply the time-split procedure of [3] in which Eq. (9) is decomposed as 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, We assume that the stress function, f S (A), and the relaxation function, f S (A), are linear functions For example, for the FENE-P constitutive model the parameters would be η S = η R = 1 and 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, 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, Step 6: Eq. (20) can be integrated explicitly, or, implicitly, 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, λ and later integrated analytically.
(a) Prior to the analytical integration, the log of the conformation tensor is diagonalised, to obtain Λ * * , R * * and the conformation tensor, Step 8: Finally,

Spatial discretization and the adaption algorithm
The open 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 [55] (see figure 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 figure 1.b, the grid would be formed by 16 leaf cells being 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 center. However, the stresses of 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.
This tree-type grid structure allows the performance of 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 in 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, This coarser field distribution, f −1 , can be upsampled (or prolongated ) to the original level, 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, 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 figure 1.c). A more detailed explanation of the adaption algorithm can be found in [56]. 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 of the prolongation and restriction operators used to decide adaption regions. In our experience, it is better to use as control adaption variables the velocity components and the volume of fraction. The values of τ p in the new refined cells are computed with a bilinear interpolation while in the coarser ones they 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 [57] 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 centers of the cells, some care must be taken to suitably model the presence of walls and symmetries in the momentum equation (16). 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 [58]. 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.

Test
We have performed various test of the numerical schemes presented in this work to verify aspect as the time integration or the correct treatment of the interaction of the viscoelastic fluid with walls and interfaces. Those tests unrelated specifically to the 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 in 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 figure 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 [27,12]. As in the previous work of [27] the dimensionless parameters were fixed to: F r = 2.26, h = 2, De = 1, Re = 5 and β = 0.1. [27] 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 Figure 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 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 (adaption A1) and ε th f = 10 −3 and ε th ux = ε th uy = 10 −2 (adaption A2). The simulation performed by [27] 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 [27] negligible difference 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 . Figure 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 [27] with the adaptation strategy A1. All three methodologies give very consistent results and are in very good agreement with the results of [27].  6). The numerical simulation of [27] 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 [51]. 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 [51]). 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 figure 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 [59], 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 [59]. 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 being the downward dimensionless velocity of the viscoelastic fluid uniform and equal to u z = −1. The rest of the variables are set to zero.
To explore the influence of the 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). Figures 5 and 6 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 of the volume fraction , c, at the wall. The effect of wall-fluid interaction is shown in Figure 5 where we plot the shape of the lamella at instant t = 0.18. The sliding lamella (red and black 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 levitating sheet is obtained by imposing the Dirichlet condition, c = 1 (green and cyan  lines). Interestingly, the elastic effects are negligible for these very low polymer concentrations as it can be observed in Figures 5 and 6. The mechanism of the splashing is unaffected by the viscoelastic character of the fluid, at least for the very small concentrations cases.Interestingly, the experiments performed by Jung et al. [60] on the splashing of droplets of solutions of polystyrene in diethyl phthalate over a highly wettable solid exhibits 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 [59], similarly to the experiments of Vega & Castrejon-Pita [51]. 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 suggest that the viscoelasticity could alter the contact line equilibrium that, in turn, affects the dynamics of the lamella. This numerical result agrees well with the experimental results of [60] 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 spreading of the droplet is unaffected by the presence of polymers. This is not the general case since the substrate will play a relevant (non simple) role in the spreading and receding stages as the dynamic contact angle will vary in the process [50,49].

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 rest of variable 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. Figure 7 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 Figure 7 also shows the 2-norm of the conformation tensor Ψ, ||Ψ|| 2 is used to visualize where the viscoelastic stresses are more intense. As can be seen in the figure 7, 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. [61] and confirmed experimentally by [62]. 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 figure 7), 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 figure  8A, 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 8C. Interestingly, the vortices roll up the viscoelastic stresses giving some sort of toroidal spring that delays the advance of the lamella (see figures 8A,B and D). 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 figure 8C. 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 figure  8D. 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. Figure D shows the details of the structure of the roll up of viscoelastic stresses around ring bubbles.

Conclusions
In this article we have shown how the time-splitting scheme proposed by Hao & Pan [3] can be used together with the classical log conformation tensor of Fattal & Kupferman [2], or the square-root conformation of Balci et al. [20], 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 [5].
The numerical scheme has been used to investigate numerically the splashing of weakly viscoelastic droplets on to solid flat substrates and pools of the same fluid, taking as reference the experimental conditions of the work of Vega & Castrejon-Pita [51]. 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. [61]. 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 some sort of toroidal spring. These toroidal springs can nucleate around trapped bubble rings similar to those reported by [61]. 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 − 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, 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, 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 figure C.10 illustrates a comparison between the analytical solution given by Eq. (C.1) with the numerical results obtained with the coarsest grid. As it can be observed, the agreement is very good and comparable to similar schemes [30], 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 figure C.11. 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) [63]. For L 2 → ∞ both FENE-CR and FENE-P coincide with Oldroyd-B. However, the plots in Fig. C.11 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  viscosity, ρ 2 and µ 2 , respectively. The interfacial surface tension is σ. Both fluids are trapped, as shown, in a planar gap of a width equal to eight times the droplet radius H = 8a, and a length approximately sixteen times the droplet radius L = 16a. Suddenly a Couette flow is imposed to both fluids, u * x (x * , y * ; t * = 0) =γ y * beingγ = 2U/H .
where the superscript * denotes dimensional variables. The rest of 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,.
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 [65] and used as a test problem by many others, see for example [27,66,67,11]. In Chinyoka et al. [65] 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. Fig. C.13A shows how this parameter evolves in our simulation (black continuous line labelled as 'Basilisk') compared with Figueiredo et al. [27] (open circles). Also, Fig. C.13B shows the position of the interface for both simulations. As expected, the agreement between both simulations is excellent.

Appendix C.3. lid cavity flow
This test deals with the movement of a viscoelastic Oldroyd-B fluid of density ρ, relaxation parameter λ, and solvent and polymeric viscosities, µ s and µ p , respectively. As shown in the insert of Fig C.14, the fluid is confined in  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 [68] 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. Figure C.14 shows the time evolution of the total dimensionless kinetic energy in the cavity, The agreement in the velocity profiles between the different methodologies, and the previous work shown in Fig. C.16, is excellent, although, this is a common result in other schemes. The agreement in the kinetic energy is also