Contact angle leads to unexpected behavior away from the wall

[Furkan Oz]

Hi,

I recently started use Basilisk in a project. We are trying to regenerate the results from this paper ( Coalescence of diffusively growing gas bubbles | Journal of Fluid Mechanics | Cambridge Core  ). However, we are having some issues with the contact angle. I started from "bubble.c" example and updated the density and viscosity ratios as well as the dimensionless numbers. Our purpose is to use initialize the bubbles from the wall, let them merge, and investigate the rising similar to paper. However, if we turn on contact angle, the behavior is unexpected and the bubble is never leaving to surface whatever we do.

In order to isolate the issue, I moved the bubbles to the center of the domain away from the walls. If I comment out the `contact_angle` lines, everything is working as expected. The moment I use `contact_angle` the bubble starts to fall instead of rising. I tried single bubble, double bubble, finer mesh, coarser mesh, different angles, different vertical direction (y and z). The behavior is always same. If I have contact angle regardless of the distance from the wall, the bubble always falls rather than rising. 

I am sorry my organization is not allowing me to attach a file but here is a copied case setup.

In this code, if I comment out contact angle boundary conditions, everything is working fine but if they are not commented, the bubbles start to fall. 

It would be great if you can show me the possible issue.

```

#include "navier-stokes/centered.h"

#include "navier-stokes/perfs.h"

#include "two-phase.h"

#include "tension.h"

/**

We can control the maximum runtime. */

#include "maxruntime.h"

#define RHOR 79.109375

#define MUR 64.26666666667

const double Ga = 2338.000161508927;

const double Bo = 0.06015151073253832;

const double MAXTIME = 110;

const double WIDTH = 10. [1];

// const double Zi = 3.5;

const double Zi = 3.5;

int LEVEL = 9;

#include "contact.h"

// #include "curvature.h"

vector h[];

h.t[bottom] = contact_angle ((90.0)*pi/180.);

h.r[bottom] = contact_angle ((90.0)*pi/180.);

u.n[bottom] = dirichlet(0.);

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

u.r[bottom] = dirichlet(0.);

/**

The main function can take two optional parameters: the maximum level

of adaptive refinement (as well as an optional maximum runtime). */

int main (int argc, char * argv[]) {

  maxruntime (&argc, argv);

  if (argc > 1)

    LEVEL = atoi (argv[1]);

 

  /**

  We set the domain geometry and initial refinement. */

 

  size (WIDTH);

  origin (-L0/2, -L0/2, 0);

  init_grid (128);

  /**

  We set the physical parameters: densities, viscosities and surface

  tension. */

 

  rho1 = 1. [0];

  rho2 = rho1/RHOR;

  mu1 = 1./Ga;

  mu2 = 1./(MUR*Ga);

  f.sigma = 1./Bo;

  f.height = h;

  /**

  We reduce the tolerance on the divergence of the flow. This is

  important to minimise mass conservation errors for these simulations

  which are very long. */

 

  TOLERANCE = 1e-4 [*];

  run();

}

/**

For the initial conditions, we first try to restore the simulation

from a previous "restart", if this fails we refine the mesh locally to

the maximum level, in a sphere of diameter 1.5 around the bubble. We

then initialise the volume fraction for a bubble initially at (0,Zi,0)

of diameter unity. */

event init (t = 0) {

  if (!restore (file = "restart")) {

    // refine (sq(x) + sq(y - Zi) + sq(z) - sq(0.75) < 0 && level < LEVEL);

    refine ((sq(x - (-0.5)) + sq(y) + sq(z - Zi) < sq(1.25) ||

           sq(x - 0.5) + sq(y) + sq(z - Zi) < sq(1.25) ||

           (fabs(x) < 3.*0.03 &&

            sq(y) + sq(z - Zi) < sq(3.*0.025))) && level < LEVEL);

    // fraction (f, sq(x) + sq(y - Zi) + sq(z) - sq(.5));

    fraction (f,

            min(min(sq(x - (-0.5)) + sq(y) + sq(z - Zi) - sq(.5),

                    sq(x - 0.5) + sq(y) + sq(z - Zi) - sq(.5)),

                max(sq(y) + sq(z - Zi) - sq(0.025),

                    sq(x) - sq(0.03))));

  }

}

/**

We add the acceleration of gravity (unity) in the downward (-y/-z)

direction. */

event acceleration (i++) {

  face vector av = a;

  foreach_face(z)

    av.z[] -= 1.;

}

/**

We adapt the mesh by controlling the error on the volume fraction and

velocity field. */

event adapt (i++) {

  double uemax = 1e-2;

  adapt_wavelet ({f,u}, {0.01,uemax,uemax,uemax}, LEVEL, 5);

}

/**

## Outputs

Every ten timesteps, we output the time, volume, position, and

velocity of the bubble. */

event logfile (i += 10) {

  double xb = 0., yb = 0., zb = 0., sb = 0.;

  double vbx = 0., vby = 0., vbz = 0.;

  foreach(reduction(+:xb) reduction(+:yb) reduction(+:zb)

    reduction(+:vbx) reduction(+:vby) reduction(+:vbz)

    reduction(+:sb)) {

    double dv = (1. - f[])*dv();

    xb += x*dv;

    yb += y*dv;

    zb += z*dv;

    vbx += u.x[]*dv;

    vby += u.y[]*dv;

    vbz += u.z[]*dv;

    sb += dv;

  }

  fprintf (stderr,

     "%.8f %.8f %.8f %.8f %.8f %.8f %.8f %.8f\n",

     t, sb,

     xb/sb, yb/sb, zb/sb,

     vbx/sb, vby/sb, vbz/sb);

  fflush (stderr);

}

event logfile (i += 1; t <= 1) {

  double ub = 0., ui = 0.;

  foreach (reduction(max:ub) reduction(max:ui)) {

    double um = sqrt(sq(u.x[]) + sq(u.y[]) + sq(u.z[]));

    if (z < 0.5) ub = max(ub, um);        

    if (fabs(z - Zi) < 1.) ui = max(ui, um);

  }

  fprintf(stderr, "%g %g %g\n", t, ub, ui);

}

#include "output_vtu_foreach.h"

event paraview (i += 250)

{

  char subname[80], fname[80];

  sprintf(subname, "bubble-%8.8d", (int) i);

#if _MPI

  sprintf(fname, "%s_n%3.3d.vtu", subname, pid());

#else

  sprintf(fname, "%s.vtu", subname);

#endif

  FILE * fp = fopen(fname, "w");

  output_vtu_bin_foreach((scalar *){f, p}, (vector *){u}, fp, false);

  fclose(fp);

#if _MPI

  if (pid() == 0) {

    sprintf(fname, "%s.pvtu", subname);

    fp = fopen(fname, "w");

    output_pvtu_bin((scalar *){f, p}, (vector *){u}, fp, subname);

    fclose(fp);

  }

#endif

}

event end (t = 0.5)

{

  if (pid() == 0)

    fprintf (stderr, "Finished at t=%g s\n", t);

} ```