# Impact of a viscoelastic drop on a solid

We solve the axisymmetric, incompressible, variable-density, Navier–Stokes equations with two phases and use the log-conformation method to include viscoelastic stresses. The curvature module is used to compute interface properties.

``````#include "axi.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "log-conform.h"
#include "curvature.h"
``````

The density and viscosity ratios are 1000. The Reynolds number based on the droplet diameter, velocity and viscosity is 5 and the Froude number is 2.26.

``````#define RHO_r 0.001
#define MU_r 0.001
#define RE 5.
#define FR 2.26
#define LEVEL 7``````

The dimensionless viscoelastic properties used are the ratio of the solvent to the total viscoelastic viscosity (polymeric plus solvent) and the Weissenberg number.

``````#define BETA 0.1
#define WI 1.0

scalar lambdav[], mupv[];``````

The drop comes from the right. We allow the fluid to get through that boundary.

``````u.n[right] = neumann(0);
p[right]   = dirichlet(0);``````

The wall is at the left side. We apply a no-slip boundary condition and a non-wetting condition for the VOF tracer.

``````u.t[left] = dirichlet(0);
tau_qq[left] = dirichlet(0);
f[left]   = 0.;

int main() {``````

The domain spans $\left[0:2.6\right]×\left[0:2.6\right]$.

``````  size (2.6);
init_grid (1 << LEVEL);``````

The densities and viscosities are defined by the parameters above.

``````  rho1 = 1.;
rho2 = RHO_r;
mu1 = BETA/RE;
mu2 = MU_r/RE;``````

The viscoelastic fields will be set below.

``````  mup = mupv;
λ = lambdav;``````

We set a maximum timestep. This is necessary for proper temporal resolution of the viscoelastic stresses.

``````  DT = 2e-3;
run();
}

event init (t = 0) {``````

At a wall of normal $\mathbf{\text{n}}$ the component of the viscoelastic stress tensor $ta{u}_{{p}_{nn}}$ is zero. Since the left boundary is a wall, we set $ta{u}_{{p}_{xx}}$ equal to zero at that boundary.

``````  scalar s = tau_p.x.x;
s[left] = dirichlet(0.);``````

The drop is centered on (2,0) and has a radius of 0.5.

``  fraction (f, - sq(x - 2.) - sq(y) + sq(0.5));``

The initial velocity of the droplet is -1.

``````  foreach()
u.x[] = - f[];
}``````

We add the acceleration of gravity.

``````event acceleration (i++) {
face vector av = a;
foreach_face(x)
av.x[] -= 1./sq(FR);
}``````

We update the viscoelastic properties. Only the droplet is viscoelastic.

``````event properties (i++) {
foreach() {
mupv[] = (1. - BETA)*clamp(f[],0,1)/RE;
lambdav[] = WI*clamp(f[],0,1);
}
boundary ({mupv, lambdav});
}``````

We adapt the solution at every timestep based on the interface and velocity errors.

``````#if TREE
adapt_wavelet ({f, u.x, u.y}, (double[]){1e-2, 5e-3, 5e-3},
maxlevel = LEVEL, minlevel = LEVEL - 2);
}
#endif``````

We track the spreading diameter of the droplet.

``````event logfile (i += 20; t <= 5) {
scalar pos[];
position (f, pos, {0,1});
fprintf (stderr, "%g %g\n", t, 2.*statsf(pos).max);
}``````

We generate a movie of the interface shape.

``````#include "view.h"

event viewing (i += 10) {
view (width = 400, height = 400, fov = 20, ty = -0.5,
quat = {0, 0, -0.707, 0.707});

clear();
draw_vof ("f", lw = 2);
squares ("u.x", linear = true);
box (notics = true);
mirror ({0,1}) {
draw_vof ("f", lw = 2);
squares ("u.y", linear = true);
box (notics = true);
}
save ("movie.mp4");``````

We can optionally visualise the results while we run.

``````#if 0
static FILE * fp = popen ("bppm","w");
save (fp = fp);
#endif
}``````

## Results

Animation of the interface shape. The color field on the right-hand-side (resp. l.h.s.) is the radial (resp. axial) velocity component.

The time evolution of the maximum diameter is close to that reported by Figueiredo et al. (2016).

## References

 [figueiredo2016] RA Figueiredo, CM Oishi, AM Afonso, IVM Tasso, and JA Cuminato. A two-phase solver for complex fluids: Studies of the Weissenberg effect. International Journal of Multiphase Flow, 84:98-115, 2016.