Viscoelastic 2D drop in a Couette Newtonian shear flow

This problem has been used as benchmark (see for example Chinyoka et al. (2005)). Equations are usually made dimensionless with the outer density ρ2, the droplet radius a and the shear rate γ̇. The following dimensionless parameters appear We=ρ2a3γ̇2σ,Re=ρ2a2γ̇μ2,μr=μ1μ2,m=ρ1ρ2,β=μsμ1,De=γ̇λ where subscript “1” and “2” stand for the droplet and matrix fluid respectively and μs is the solvent viscosity. Note that the polymeric viscosity is μp=μ1μs.

#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "log-conform.h"
#include "tension.h"

#define Ca 0.6     // Capillary number
#define Re 0.3     // Reynold number
#define We (Ca*Re) // Weber number
#define MUr 1.     // ratio of outer(matrix) to inner(drop) viscosity
#define M 1.       // ratio of outer to inner density
#define Deb 0.4    // Deborah number
#define Β 0.5   // ratio of the solvent visc. to the total viscoelastic visc.

We set a maximum level of 8 only so that the test case runs in less than 30 minutes but note that the results presented below are for MAXLEVEL = 9.

int MAXLEVEL = 8;

scalar mupd[], lam[];

The top and bottom boundary conditions are those of a Couette flow.

u.t[top] = dirichlet (y);
u.t[bottom] = dirichlet (y);

The domain is a 16x16 box which will be masked later to a bottom-left corner of the domain of coordinates (x=8,y=4). We set a maximum timestep of 0.1

int main() {
  L0 = 16 ;
  origin (-L0/2, -L0/4.);
  periodic (right);
  DT = .1;

We set the viscosities, densities, surface tension and visco-elastic parameters according to the analysis above.

  mu1 = MUr*Β/Re;
  mu2 = 1./Re;
  rho1 = M;
  rho2 = 1.;
  f.σ = 1./We;
  λ = lam;
  mup = mupd;

  init_grid (1 << MAXLEVEL);

event init (i = 0) {

We mask above y>4. The computational domain is now a 16x8 rectangle with the origin in the center.

  mask (y > 4 ? top : none);

As initial conditions we set a viscoelastic droplet of radius 1 and a linear velocity profile typical of the Couette flow.

  fraction (f, 1. - (sq(x) + sq(y)));
    u.x[] = y;
  boundary ((scalar *){u});

The event below set viscoelastic properties at each step. The dimensionless relaxation parameter is the Deborah number, De, while the dimensionless polymeric viscosity μp is μpρ2a2γ̇=μ1βμ1ρ2a2γ̇=1βReμ1μ2

event properties (i++) {
  foreach () {
    lam[] = Deb*f[];
    mupd[] = MUr*(1. - Β)*f[]/Re;
  boundary ({lam, mupd});

The mesh is adapted according to the errors on volume fraction and velocity.

event adapt (i++) {
  adapt_wavelet ({f, u}, (double[]){1e-2, 1e-3, 1e-3}, MAXLEVEL, MAXLEVEL - 2);

As outputs we plot the shape of the interface at instant t = 10 and the time evolution of the deformation parameter.

event interface (t = 10.) {
  output_facets (f);

event deformation (t += 0.1) {
  double rmax = -HUGE, rmin = HUGE ;
  foreach (reduction(max:rmax) reduction(min:rmin)) 
    if (f[] > 0 && f[] < 1) {
      coord p;
      coord n = mycs (point, f);
      double α = plane_alpha (f[], n);
      plane_area_center (n, α, &p);
      double rad  = sqrt(sq(x + Δ*p.x) + sq(y + Δ*p.y)); 
      if (rad > rmax)
	rmax = rad;
      if (rad < rmin)
	rmin = rad;
  double D = (rmax - rmin)/(rmax + rmin);
  fprintf (stderr, "%g %g %g %g\n", t, rmin, rmax, D);

We can optionally visualise the results with Basilisk View while we run.

#if 0
#include "view.h"

event viewing (i += 10) {
  static FILE * fp = popen ("bppm","w");
  view (width = 600, height = 300, fov = 10);
  draw_vof ("f");
  // squares ("u.y", linear = true);
  squares ("level");
  save (fp = fp);

We compare the results to those of Figueiredo et al. (2016).

Time evolution of the deformation

Time evolution of the deformation

Interface shape at t = 10

Interface shape at t = 10



T Chinyoka, YY Renardy, M Renardy, and DB Khismatullin. Two-dimensional study of drop deformation under simple shear for Oldroyd-B liquids. Journal of Non-Newtonian Fluid Mechanics, 130(1):45-56, 2005.


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.