sandbox/tavares/test_cases/impact-thin-fiber.c

#include "grid/multigrid3D.h"
//#include "grid/octree.h" 

#include "embed-navier.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"
#include "tavares/contact-embed.h"
#include "view.h"

#define RHOR 1000
#define MUR 50
#define T 4.
#define Rd 0.0009
#define Rf 0.000375*1.1
#define RATIO Rf/Rd
#define Df Dd*RATIO
#define Yi 0.9
#define shift 0.035*RATIO
#define yc 0.5

#if 1 // First case (the droplet is captured by the fiber)
const double Red=175;
const double Wed=1.6;
const double Frd=2;
#else
const double Red=830; // Second case (the droplet detaches from the fiber)
const double Wed=37;
const double Frd=10;
#endif

const double Dd=1.;

double theta0=15; //instead of $10^\circ$ in the reference papers

int MAXLEVEL=6;

int main() {

  size (4);
  origin (-L0/2, -L0/2, -L0/2);
  N = 1 << MAXLEVEL;
  init_grid (N);

  rho1 = 1. [0];  // 
  rho2 = rho1/RHOR;
  mu1 = 1./Red;
  mu2 = 1./(MUR*Red);
  f.sigma = 1./Wed;

 // /*
 //  We reduce the tolerance on the Poisson and viscous solvers to
 //  improve the accuracy. */
 //  TOLERANCE = 1e-4 [*];;

  const scalar c[] = theta0*pi/180.;
  contact_angle = c;
  run();
}


event init (t = 0) {

  solid (cs, fs, sq(x) + sq(y-yc) - sq(Df/2.));
  fractions_cleanup (cs, fs);

  #if ADAPT 
    refine (sq(x) + sq(y) + sq(z) < sq(1.1*Dd/2.) && level < MAXLEVEL);
  #endif

  fraction (f, - (sq(x-shift) + sq(y - (yc + Yi)) + sq(z) - sq(Dd/2.)) );

  foreach()
    u.y[] = -1*f[];
}

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

uf.n[right] = 0.;
uf.n[left] = 0.;

// u.n[embed] = dirichlet(0.);
// u.t[embed] = dirichlet(0.);
// u.r[embed] = dirichlet(0.);

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

 event logfile (i++; t<=T) {
  norm n = normf(u.y);
  fprintf(stdout, "%g %g %g %g\n", t, dt, n.avg, n.max);

 }


#if 0
event cleanfghost (i++, t<=T, last) { // clean fluid ghost cells
  foreach()
    if (cs[]==0) f[]=0.;
}
#endif


#if 1
event movie (i+=10; t <= T)
{
view (quat = {-0.058, 0.184, 0.021, 0.981},
      fov = 30, near = 0.01, far = 1000,
      tx = -0.021, ty = 0.094, tz = -3.041,
      width = 1140, height = 678, bg = {1,1,1});
  //cells (n = {0,1,0});
  draw_vof (c = "cs", s = "fs", filled = -1, fc = {0.83,0.83,0.83});
  draw_vof (c = "f", fc = {0.647,0.114,0.176}, lw = 2);
  save ("movie.mp4");
}
#endif

#if TREE
event adapt (i++) {
fractions_cleanup(cs,fs,smin = 1.e-10);
#if 1
  scalar f1[];
  foreach()
    f1[] = f[];
  adapt_wavelet ({f1,cs}, (double[]){1e-3,1e-2}, minlevel = 3, maxlevel = MAXLEVEL);
#else
  adapt_wavelet ({f}, (double[]){1e-4}, minlevel = 3, maxlevel = MAXLEVEL);
#endif
}
#endif