sandbox/ecipriano/test/pinning.c
Pinning a Liquid Droplet
We want to test the module pinning.h, to suspend a liquid droplet at a specific point of the domain. The aim is to compute the equilibrium contact angle and to compare it with the analytical solution:
\displaystyle \theta_{2D} = \arccos \left( \dfrac{\rho_l g \pi r_d^2}{2 \sigma}\right)
\displaystyle \theta_{AXI} = \arccos \left( \dfrac{4\rho_l g r_d^3}{3 \sigma d_f}\right)
given by the balance between gravity and surface tension force.
#if AXI
# include "axi.h"
#endif
#include "navier-stokes/centered.h"
#include "pinning.h"
#define FILTERED
#include "two-phase.h"
#include "tension.h"
#include "reduced.h"
#include "view.h"
Boundary Conditions
We set the velocity to zero on the bottom
of the domain, in order to mimic the effect of a solid fiber.
u.n[bottom] = dirichlet (0.);
u.t[bottom] = dirichlet (0.);
p[bottom] = neumann (0.);
uf.n[bottom] = 0.;
uf.t[bottom] = 0.;
Simulation Setup
We initialize a volume fraction fn
which is used to evaluate the convergence of the numerical simulation, based on the small variation of volume fraction between two consecutive time steps.
#define CONVERGENCE_TOLERANCE 1.e-10
scalar fn[];
We set the diameter of the droplet and the maximum level of refinement.
int maxlevel = 7;
double D0 = 1.e-3 [*];
int main (void) {
We set the density and viscosity to water-air properties, in order to ensure density and viscosity ratios which are representative of realistic situations.
rho1 = 1000., rho2 = 1.;
mu1 = 1.e-3, mu2 = 1.e-5;
We set the length of the domain, and we move the origin according to the fiber diameter along y (specific to AXI).
L0 = 4.*D0;
X0 = -0.5*L0;
#if AXI
Y0 = 0.08*D0;
#endif
We set the gravity contribution along the x direction according to Basilisk convention for AXI simulations.
G.x = -9.81;
Data for the pinning
model: we set the pinning point to the initial radius of the droplet, considering the fiber diameter in AXI coordinates.
pinning.ap = sqrt (sq (0.5*D0) - sq (Y0));
pinning.ac = 0.;
We run the simulation at different values of surface tension. Shallower surface tension values would require the implementation of the normal components of the height function as explained in contact.h.
#if AXI
for (f.sigma = 0.016; f.sigma <= 0.06; f.sigma += 0.004)
#else
for (f.sigma = 0.006; f.sigma <= 0.06; f.sigma += 0.004)
#endif
{
init_grid (1 << maxlevel);
run();
}
}
We initialize half liquid droplet on the bottom boundary.
#define circle(x,y,R)(sq(R) - sq(x) - sq(y))
event init (t = 0) {
fraction (f, circle (x, y, 0.5*D0));
foreach()
fn[] = f[];
}
We make sure that the interface is never unrefined.
event adapt (i++)
{
scalar fa[];
foreach() {
if (interfacial (point, f))
fa[] = rand();
else
fa[] = 0.;
}
adapt_wavelet ({fa,u.x,u.y}, (double[]){1.e-3,1.e-3,1.e-3}, maxlevel);
}
Post-Processing
The following lines of code are for post-processing purposes.
Maximum Velocity
We output the maximum velocity to verify that the velocity relaxes toward a null value.
event logger (i += 100)
{
scalar un[];
foreach()
un[] = norm(u);
fprintf (stdout, "%g %g\n", t, normf(un).max);
}
Contact Angle and Droplet Shape
We write a function that calculates and write the numerical contact angle from the height–functions, and we write the droplet shape using the facets.
void write_theta (void) {
double ca = 0.;
foreach_point (pinning.ap, Y0)
ca = atan (1./(h.x[0,-1] - h.x[]));
fprintf (stdout, "\n\n");
fflush (stdout);
fprintf (stderr, "%g %g\n", f.sigma, fabs (ca*180./pi));
fflush (stderr);
}
void write_facets (void) {
char name[80];
sprintf (name, "facets-%.3f", f.sigma);
FILE * fp = fopen (name, "w");
output_facets (f, fp);
fclose (fp);
}
We stop the simulation if the variation of volume fraction between two consecutive time steps is smaller than 1.e-10, similarly to what is done in spurious.c.
event stop (i++)
{
double dc = change (f, fn);
if (i > 1 && dc < CONVERGENCE_TOLERANCE) {
write_theta();
write_facets();
return 1;
}
}
At equilibrium (t = 10 seems sufficient), we compute the numerical contact angle.
event end (t = 2.)
{
write_theta();
write_facets();
}
Results
We compare the steady state shape of the droplet at different surface tension values.
reset
set term push
set term @SVG size 640,180
set size ratio -1
set xrange[-1.1e-3:1.1e-3]
set yrange[0:0.5e-3]
unset xtics
unset ytics
unset border
plot 0 lt -1 notitle, \
"facets-0.006" w l lw 2 t "sigma = 0.006", \
"facets-0.010" w l lw 2 t "sigma = 0.010", \
"facets-0.014" w l lw 2 t "sigma = 0.014", \
"facets-0.030" w l lw 2 t "sigma = 0.030", \
"facets-0.058" w l lw 2 t "sigma = 0.058"
set term pop
reset
set term push
set term @SVG size 640,180
set size ratio -1
set xrange[-1.1e-3:1.1e-3]
set yrange[0:0.5e-3]
unset xtics
unset ytics
unset border
plot 0.08*1.e-3 lt -1 notitle, \
"../pinning-axi/facets-0.016" w l lw 2 t "sigma = 0.016", \
"../pinning-axi/facets-0.020" w l lw 2 t "sigma = 0.020", \
"../pinning-axi/facets-0.024" w l lw 2 t "sigma = 0.024", \
"../pinning-axi/facets-0.040" w l lw 2 t "sigma = 0.040", \
"../pinning-axi/facets-0.056" w l lw 2 t "sigma = 0.056"
set term pop
We plot the comparison between the equilibrium contact angle and the one obtained from the numerical simulation.
reset
set yrange[20:100]
set xrange[0:0.06]
set xlabel "Gravity"
set ylabel "Contact Angle [deg]"
set grid
f(x) = acos(pi*(0.5e-3)**2*1000.*9.81/2/x)*180/pi # 2D
plot f(x) w l lw 2 t "Theoretical", \
"log" u 1:2 w p pt 8 ps 1.5 t "Numerical"
reset
set yrange[0:100]
set xrange[0:0.06]
set xlabel "Gravity"
set ylabel "Contact Angle [deg]"
set grid
df = 2.*0.08*1e-3
f(x) = acos(4.*1000*9.81*(0.5e-3**3)/(3*x*df))*180/pi # 3D
plot f(x) w l lw 2 t "Theoretical", \
"../pinning-axi/log" u 1:2 w p pt 8 ps 1.5 t "Numerical"
We plot the trend of the maximum velocity which relaxes toward a null value during the simulation time.
reset
set xrange[0:2]
set xlabel "time [s]"
set ylabel "Maximum velocity norm [m/s]"
set grid
plot "out" i 0 w l t "sigma = 0.006", \
"out" i 4 w l t "sigma = 0.022", \
"out" i 6 w l t "sigma = 0.03" , \
"out" i 13 w l t "sigma = 0.058"
reset
set xrange[0:2]
set xlabel "time [s]"
set ylabel "Maximum velocity norm [m/s]"
set grid
plot "../pinning-axi/out" i 0 w l t "sigma = 0.016", \
"../pinning-axi/out" i 4 w l t "sigma = 0.032", \
"../pinning-axi/out" i 7 w l t "sigma = 0.044" , \
"../pinning-axi/out" i 10 w l t "sigma = 0.056"