sandbox/Antoonvh/axibounce.c

Astronaut Scott Kelly played ping pong in space. Video courtesy of NASA Johnson (\leftarrow click to watch the full movie)

Astronaut Scott Kelly played ping pong in space. Video courtesy of NASA Johnson ( click to watch the full movie)

A Bouncing Axisymetric Droplet

On this page a quasi three dimensional (3D) simulation of the bouncing droplet is presented. The set-up is very similar to that of the planar (2D) version here, but now in the 3D-ish, axisymmetric limit.

#include "axi.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"

The case set-up

Our “vertical” axis, if that concept is suitable to use inside a space station, is now the swapped with the former horizontal axis. Hence, what used to be associated with the bottom boundary in the 2D example is now called left. Furthermore, the radial coordinate (r) is associated with the y-direction in Basilisk’s axisymmetric formulation.

f[left]=0.;
u.t[left]=dirichlet(0.);
f[right]=0.;
u.t[right]=dirichlet(0.);

The physical system consists of a liquid droplet with a radius R that travels trough the air towards the hydrophobic wall with a velocity U. Gauged from the movie, we approximate R2cm and U3cm/s.Togetherwiththefluidpropertiesoftheairandwater(a,w,a,w, $) we can identify some dimensionless groups that describe the system.

Re=ρaURμa30,

We=ρwU2Rσ0.2,

Π1=μwμa100,

Π1=ρwρa1000100,

D=32; Axisymmetric,

where D stands here for the number of spatial dimensions of the system and the ‘’ symbols indicate dimensionless groups where we make a consession with respect to the reality onboard the space station to keep the numerical experiment feasable on a single-core system.

We set the model parameters such that we obtain the approximated ratios, aproximately. It is done in such a way that we have a normalized velocity scale (U), length scale (R) and gas phase density (ρa) scale. Furthermore, we set the maximum grid resolution to correspond to a 2562 equidistant grid. For some minimally desired consistency, we take special care such that the radial coordinate does not have any negative values.

double R = 1.;
double U = 1.;
int maxlevel = 8;
double temp = 50.;

int main(){
  mu2=1./30.; //gas phase
  mu1=100./30.; //liquid phase
  rho2=1.; //gas phase
  rho1=100;//liquid phase
  f.σ=500.;
  init_grid(1<<7);
  L0=10;
  X0=-L0/2;
  Y0=0.;
  run();
}

initialization

After we have refined the grid with a ring of high resolution mesh, the droplet is initialized and it is targeted at the left wall.

event init(i=0){
  refine(sq(x)+sq(y)<sq(R+0.25) && sq(x)+sq(y)>sq(R-0.25) && level<maxlevel);
  fraction(f,sq(R)-sq(x)-sq(y));
  foreach()
    u.x[]=-f[]*U;
}

## Adaptation Since the advective interface tracking and resolving for the surface tencile waves only requires a high resolution mesh at the locations of the interface, it seems smart to use grid adaptivity, we will check it this was indeed the case.

event adapt(i++)
  adapt_wavelet((scalar *){u,f},(double []){0.02,0.02,0.001},maxlevel);

Output

The general dynamics are visualized in a movie that shows the rebound of the droplet.

event view(t=0.05;t+=0.1;t<=temp){
  static FILE * fp5 = popen("gfsview-batch2D axibounce.justphase.gfv | ppm2mp4 axiphase.mp4","w"); 
  output_gfs(fp5);
  fprintf (fp5, "Save stdout { format = PPM width = 800 height = 450}\n");
}

The simulation does display the ping pong dynamics;

The energy partitioning between kinetic and the surface free energy is again calculated. The formulation as was used in the planar simulation is modified to correct for the fact that the y coordinate now represents the radial coordinate.

void find_facets(Point point,scalar f,double xy[4]){
  coord n;
  n = mycs (point, f);
  double α = plane_alpha (f[], n);
  coord segment[2];
  if (facets (n, α, segment) == 2){
    xy[0] = x + segment[0].x*Δ;
    xy[1] = y + segment[0].y*Δ;
    xy[2] = x + segment[1].x*Δ;
    xy[3] = y + segment[1].y*Δ;
  }else{
    printf("Warning:\nCould not find facets; expect unexpected behaviour.\n");
  }
}

event energy(i+=5){
  static FILE * fpe = fopen("axienergy","w");
  double e=0;
  double a=0;
  double xyf[4];
  foreach(reduction(+:e) reduction(+:a)){
    e+=0.5*(sq(Δ)*M_PI*2*y)*(sq(u.x[])+sq(u.y[])) * ((rho1*f[]) + (rho2*(1-f[])));
    if (f[]>0.00001 && f[] < 0.9999){
      find_facets(point,f,xyf);
      a+=2*M_PI*((xyf[1]+xyf[3])/2.)*pow(sq(xyf[0]-xyf[2])+sq(xyf[1]-xyf[3]),0.5);
    }
  }
  fprintf(fpe,"%g\t%d\t%g\t%g\t%g\t%g\n",t,i,e,(a-(4*M_PI*sq(R)))*f.σ,a,((a-(4*M_PI*sq(R)))*f.σ) + e);
  printf("%g\t%d\t%g\t%g\t%g\t%g\n",t,i,e,(a-(4*M_PI*sq(R)))*f.σ,a,((a-(4*M_PI*sq(R)))*f.σ) + e);
}

We can view the results;

Overall, things seem familiar and consistent to what was learned from the planar case. Differences may be spotted when doing a more quantitative analysis of the energy partioning.

The next step

In the movie, the droplet does not really appear to be axisymetric. It is generally quite challenging to make-up an a priori argument on how to dynamics would alter when the physical system is allotted with an additional spatial dimension. E.g. after having studied 2D turbulence (a la Kaichmann), would you be able to expect the main flow-behavioural characteristics of 3D turbulence (a la Kolmogorov)? Well, Not me! Therefore, it is required to break the axisymetry and study the flow in a fully-fleched 3D simulation. This study should also adresses the influence of the consession that was made with regards to the ratio of the phases’ densities.