sandbox/popinet/dewetting.c
Dewetting
Implemented from Mahady, Afkhami, Kondic, Phys. Fluids 28, 062002 (2016).
Several things need to be improved: adaptivity, more robust height estimation (as done by VariablePosition).
#include "grid/multigrid.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"
double R = 0.4, hs = 0.01, thetai = pi/2., Oh = 1./sqrt(10.);
int main() {
N = 128;
f.sigma = 1;
mu1 = mu2 = sq(Oh)*f.sigma*2.*R;
u.t[bottom] = dirichlet(0);
f[bottom] = 1.;
run();
}
event init (t = 0) {
#if 0
fraction (f, 5.*hs*(1. + 0.01*sin(4.*M_PI*x)) - y);
#elif 0
fraction (f, max(sq(R) - sq(x) - sq(y + R*cos(thetai) - hs),
hs - y));
#else
fraction (f, 5.*hs*(1. + 0.01*noise()) - y);
#endif
}
event acceleration (i++)
{
scalar hy[];
foreach()
f[] = clamp (f[], 0, 1);
position (f, hy, {0,1});
#if TREE
f.prolongation = p.prolongation;
f.dirty = true;
#endif
// Eq. 10 of Mahady et al, 2016
double theta_eq = pi/2., m = 3, n = 2;
double xi = (1. - cos(theta_eq))/hs*(m - 1.)*(n - 1.)/(m - n);
face vector st = a;
foreach_face()
if (f[] != f[-1]) {
double h =
(hy[] < nodata && hy[-1] < nodata) ?
(hy[] + hy[-1])/2. :
hy[] < nodata ? hy[] :
hy[-1] < nodata ? hy[-1] :
nodata;
if (h != nodata)
st.x[] -= alpha.x[]/fm.x[]*f.sigma*xi*(pow(hs/h, m) - pow(hs/h, n))
*(f[] - f[-1])/Delta;
}
#if TREE
f.prolongation = fraction_refine;
f.dirty = true;
#endif
}
#if 0
event gfsview (i += 10) {
static FILE * fp = popen ("gfsview2D -s dewetting.gfv", "w");
output_gfs (fp, t = t);
}
#endif
event logfile (i++) {
fprintf (stderr, "%g %g\n", t, statsf (u.x).max);
}
event interface (i += 200; t <= 4) {
output_facets (f);
}
Results
plot 'out' w l t '', 0.01