/**
![Water droplets may rest on a hydrophobic surface. Image via [the university of Wisconsin](https://www.chem.wisc.edu/content/hydrophobic-interactions).](https://www.chem.wisc.edu/deptfiles/Water.png)
# Equilibrium shapes of a droplet on a hydrophobic surface
An atmosphere with a density $\rho_2$ hosts a denser, watery droplet with a
density $\rho_1 = \rho_2 + \Delta\rho$. The droplet is affected by the
acceleration due to gravity ($g$) and rests on a perfect hydrophobic
surface. A characteristic lengthscale of the droplet ($R$) may be
introduced by considering it's volume; $V=4/3 \pi R^3$. We assume that
there exists a single equilibrium state of the system that is governed
by the so-called Bond number:
$$Bo = \frac{\Delta \rho g R^2 }{\sigma}$$
On this page we aim to find the shape of the droplet for a range of
$Bo$ values. Because of symmetry reasons, we choose to use the
axi-symmetric solver to save some computational effort.
*/
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"
#include "axi.h"
#include "reduced.h"
#include "view.h"
int maxlevel = 8;
double R = 1.;
f[left] = 0.;
u.t[left] = dirichlet (0.);
int sw = 1;
char fname[99];
int it = 0;
double AT = 0.5*10E-4;
/**
## set-up
We choose to use normalized values for all parameters that define the
Bond number except for the value of the surface tension
($\sigma$). This parameter is used to control the aforementioned
dimensionless ratio $Bo$.
*/
int main(){
G.x = -1;
mu1 = 1.;
mu2 = 1.;
rho1 = 2.;
rho2 = 1.;
L0 = 4*R;
init_grid (64);
run();
}
/**
For $Bo \rightarrow 0$, a spherical droplet shape is expected. We
therefore start with an small value for the Bond number and initialize
the droplet as a sphere.
*/
event init (t=0){
TOLERANCE = 10E-4;
f.sigma = 32.;
refine (sq(x-R) + sq(y) < 1.1*R && level < maxlevel);
fraction (f, R - sq(x-R) - sq(y));
}
/**
## Finding equilibrium solutions
Due to the effect of gravity, the droplet will deform. Deformation is
characterized by velocities and we assume that equilibrium solutions
have no kinetic energy in the system. We apply an arbritrary threshold
(`AT`) that sets a lower limit for the kinetic energy. When the
kinetic energy is lower than this value, the solution is considered to
be in equilibrium. At this point we output the shape of the bubble to
a file. Next, the value of the Bond number is increased by a factor of
two and a new equilibrium is calculated, starting from the
previous solution. Since this is not really a physically consistent set-up we also monitor the volume of the droplet to diagnose possible changes in the $R$ parameter.
*/
double e = 0.;
event check_equi (i+=25;t<100){
double en=0.;
foreach(reduction(+:en)){
foreach_dimension()
en+=sq(u.x[])*sq(Delta);
}
if (sw == 0){
it++;
if (it > 10)
sw=1;
}
if (en