src/test/marangoni.c

    Marangoni-induced translation due to a temperature gradient

    This reproduces the test case in section 3.4 of Al Saud et al., 2018 which should be consulted for more details.

    Final velocity field, interface, surface tension gradient and adaptive mesh

    Final velocity field, interface, surface tension gradient and adaptive mesh

    We reproduce below Figure 13.a of Al Saud et al., 2018, but for a different capillary number (also note the different vertical scale).

    set xlabel 't*'
set ylabel 'u_{drop}/U_{drop}'
set yrange [0.8:]
set key bottom right
plot 'out' i 2 u 5:6 w l t '{/Symbol D} = 1/32 R', \
     'out' i 1 u 5:6 w l t '{/Symbol D} = 1/16 R', \
     'out' i 0 u 5:6 w l t '{/Symbol D} = 1/8 R'
    Ratio of the translation velocity and theoretical velocity for different resolutions. (script)

    Ratio of the translation velocity and theoretical velocity for different resolutions. (script)

    The rate of convergence and the level of error is also good compared to Table 4 of Al Saud et al., 2018.

    reset
set xlabel 'Number of grid points per radius'
set ylabel 'Relative error'
set logscale
set xtics 8,2,32
set xrange [6:40]
set grid
plot 'log' u 1:(($3 - $2)/$3) t '', 2.3*x**-2 t 'x^{-2}'
    Rate of convergence toward the theoretical terminal velocity (script)

    Rate of convergence toward the theoretical terminal velocity (script)

    #include "axi.h"
#include "navier-stokes/centered.h"
#include "two-phase-clsvof.h"
#include "integral.h"
#include "view.h"

int LEVEL = 8;

    See section 3.4 of Al Saud et al., 2018. Note that we use a capillary number Ca 10 times larger than in Al Saud et al. to make the computation faster, but the results are good also for \text{Ca} = 0.066.

    const double R = 1. [1], NablaT = 1., Mu = 1., Rho = 1. [0];
const double Re = 0.066, Ca = 0.66;
const double Gamma_T = Re*sq(Mu)/(Rho*sq(R)*NablaT);
const double Gamma_0 = (Gamma_T*R*NablaT)/Ca;
const double t0 = Mu/(Gamma_T*NablaT);
const double Cdrop = 1., Cbulk = 1.;
double U_drop;

    We need a field to store the variable surface tension coefficient.

    scalar sigmaf[];

int main()
{
  size (16*R);
  origin (- L0/2.);
  rho1 = rho2 = Rho;
  mu1 = mu2 = Mu;
  d.sigmaf = sigmaf;

    We reduce the tolerance on the Poisson and viscous solvers to improve the accuracy.

      TOLERANCE = 1e-4 [*];
  
  U_drop = - 2./((2. + 3.*mu2/mu1)*(2. + Cdrop/Cbulk))*Gamma_T*R*NablaT/mu1;

  for (LEVEL = 7; LEVEL <= 9; LEVEL++) {
    N = 1 << LEVEL;
    run();
  }
}

    We initialize the signed distance d and the surface tension gradient.

    event init (t = 0)
{
  foreach() {
    d[] = sqrt (sq(x) + sq(y)) - R;
    sigmaf[] = Gamma_0 + Gamma_T*NablaT*x;
  }
}

    A utility function to check the convergence of the multigrid solvers.

    void mg_print (mgstats mg)
{
  if (mg.i > 0 && mg.resa > 0.)
    fprintf (stdout, "%d %g %g %g %d ", mg.i, mg.resb, mg.resa,
	     mg.resb > 0 ? exp (log (mg.resb/mg.resa)/mg.i) : 0.,
	     mg.nrelax);
}

    We log the position of the center of mass of the bubble, its velocity and volume as well as convergence statistics for the multigrid solvers.

    double u_drop = 0.;

event logfile (i += 10)
{
  double xb = 0., vb = 0., sb = 0.;
  static double xb0 = 0., previous = 0.;
  if (t == 0.)
    previous = 0.;
  foreach(reduction(+:xb) reduction(+:vb) reduction(+:sb)) {
    double dv = (1. - f[])*dv();
    vb += u.x[]*dv;
    xb += x*dv;
    sb += dv;
  }
  static double sb0 = 0.;
  if (i == 0) {
    sb0 = sb;
    fprintf (stdout, "\nt dsb xb vb/U_drop ta u_drop/U_drop dt perf.t perf.speed\n");
  }
  u_drop = t > previous ? (xb/sb - xb0)/(t - previous) : 0.;
  fprintf (stdout, "%g %g %g %g %g %g %g %g %g ",
	   t/t0, (sb - sb0)/sb0, xb/sb, vb/sb/U_drop,
	   (t + previous)/2./t0, u_drop/U_drop,
	   dt, perf.t, perf.speed);
  mg_print (mgp);
  mg_print (mgpf);
  mg_print (mgu);
  fputc ('\n', stdout);
  xb0 = xb/sb, previous = t;
  fflush (stdout);
}

event graphics (t = 3.*t0)
{
  double U = - Gamma_T*R*NablaT/Mu;
  fprintf (stderr, "%d %g %g\n", N/16, u_drop/U, U_drop/U);
  if (LEVEL == 8) {
    view (fov = 30, near = 0.01, far = 1000,
	  tx = 0.009, ty = -0.076, tz = -0.291,
	  width = 1239, height = 575);
    draw_vof (c = "f", filled = - 1, fc = {1,1,1});
    draw_vof (c = "f", lw = 2);
    squares (color = "sigmaf", spread = 0.8, linear = true);
    vectors (u = "u", scale = 1);
    cells ();
    save ("fields.png");
  }
}

#if TREE
event adapt (i++) {
  adapt_wavelet ({f,u}, (double[]){1e-2, 1e-5, 1e-5}, LEVEL);
}
#endif

    References

    [alsaud2018]

    Moataz O. Abu-Al-Saud, Stéphane Popinet, and Hamdi A. Tchelepi. A conservative and well-balanced surface tension model. Journal of Computational Physics, 371:896–913, February 2018. [ DOI | http | .pdf ]