sandbox/Antoonvh/coffee.c

Some people prefer to mix some milk in their coffee. Image via Antony Hall.

Mixing milk with Coffee

A liquid layer with depth h_0 and density \rho_1 is hosted in a less dense atmosphere (\rho_2) and kept in a mug-like geometry with radius R by the acceleration due to gravity (g). It is brought into a solid-body rotation with an angular velocity v_{\theta} = \omega. At some point the mug becomes stationary and due the the viscous boundary layers at the no slip walls a secondary (Ekman) circulation will be generated. This warrants to introduce both fluid’s viscosities (\mu_1 and \mu_2), giving rise to the following dimensional groups, with typical values for coffee mixing:

\displaystyle \frac{h_0}{R}= 0.5

\displaystyle \frac{\rho_1}{\rho_2} = 1000

\displaystyle \frac{\mu_1}{\mu_2} = 50

\displaystyle \frac{\omega^2 R}{g} = 0.2

\displaystyle \frac{\rho_1 \omega R^2}{\nu_1} = 15000

Using the approximate properties of water and air, R=5\mathrm{cm}, \omega = 2\pi\mathrm{Rad/s}, g = 9.81 \mathrm{m}^2/\mathrm{s} and a near-empty cup of coffee h_0 = 2.5\mathrm{cm}. For the dimensionless set up we choose g = 1, R = 1 and set \rho_1 = 1, then: h_0 = 0.5R, \omega = 0.45\sqrt{g/R}, \nu_1 = 3\times 10^{-5}\rho_1 \omega R^2, \nu_2 = 6\times 10^{-7} \rho_1 \omega R^2 and \rho_2 = 1/1000 \rho_1. Noting that we do not consider the surface tension associated with the liquid-atmosphere interface (\sigma) an important parameter for the dynamics (Bo \rightarrow \infty).

Numerical Set up

We assume the swirling flow to be axisymmetric, this reduces the computational costs considerably and makes it easy to implement the circular mug boundaries.

#include "grid/multigrid.h" // I do not understand axi + adaptivity
#include "axi.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "reduced.h"
#include "view.h"
#include "particles.h"

Swirl

By default the axisymmetric flow solver only considers u_z and u_r, we need to extend it with the angular velocity (in Rad/s) u_{\theta}. We advect and diffuse the angular velocity field (ut), naively assuming this works well with the axi-symmetric metric. Furtheremore, we take into account the centrifugal force due to the non-zero swirl.

#include "tracer.h"
#include "diffusion.h"
scalar ut[];
scalar * tracers = {ut}; 
face vector av[];

event acceleration(i++){
  coord f = {0, 1};
  foreach_face(y){
    if (y > Delta/2.)
      av.y[] += f.y*y*sq((ut[0, -1] + ut[])/2.);
    else
      av.y[] += 0.;
  }
  DT = (min(DT*1.01, 0.1));
}

event tracer_diffusion(i++)
  diffusion(ut, dt, mu);

#define BC (0.45*(t < tf))
double tf = 1;
uf.n[top] = 0;// This solved a prominent mass-conservation issue
ut[top] = dirichlet(BC);
ut[left] = dirichlet(BC);

u.t[left] = dirichlet(0.);
u.t[top] = dirichlet(0.);
f[bottom] = neumann(0.);

We initialize and set some more properties before we start the run().

int main(){
  Y0 = 0.;
  init_grid(128);
  a = av;
  mu1 = 3E-5;
  mu2 = 6E-7;
  G.x = -1;
  G.y = 0.;
  rho1 = 1.;
  rho2 = 1./1000.;
  run();
}

Mixing of milk

In order to check if the aforementioned secondary circularion is an efficient mixing mechanism we add tracer particles to the flow.

event init( t = 0){
  init_particles_2D_square_grid(50, 0.1, 0.10, 0.1); 
  int j = 0;
  while (j < n_part)
    loc[j++].z = 0; // This is not axisymmetric
  fraction(f, (0.5 + 0.5*sq(BC)*sq(y)) - x);
  TOLERANCE = 10E-5;
  DT = 0.001;
  foreach()
    ut[] = BC;
  boundary(all);
}

Apart from the particle advection in the r and z direction, the \theta direction needs to be threated as well.

event advance_particles(i++){
  int j = 0;
  while (j < n_part){
    loc[j].z += dt*interpolate(ut, loc[j].x, loc[j].y);
    j++;
  }
}

Generating the visialization

If we wish to view the paths of the tracers, some extra care is required not to let the 2D Basilisk view functionalities limit our artistic freedom.

// This is copied from draw.h
static void glvertex3d (bview * view, double x, double y, double z) { 
  if (view->map) {
    coord p = {x, y, z};
    view->map (&p);
    glVertex3d (p.x, p.y, p.z);
  }
  else
    glVertex3d (x, y, z);
}

struct _scatter {
  coord * loc;    // Coordinates
  float s, pc[3], coefs[3]; // point size, colour and distance attenuation coefs.
};

void glPointParameterfv (GLenum pname, const GLfloat * params);

trace
void scatter3Din2D(struct _scatter p){// Modified from scatter.h
  bview * view = draw();
  if (!p.coefs[0]){ // A guess:
    p.coefs[0] = 0.01;
    p.coefs[1] = 0.2;
    p.coefs[2] = 0.5;
  }
  glPointParameterfv(GL_POINT_DISTANCE_ATTENUATION, p.coefs);
  if (p.pc)
    glColor3f(p.pc[0], p.pc[1], p.pc[2]);
  if (!p.s)
    p.s = 20;
  glPointSize(p.s);
  glBegin (GL_POINTS);
  for (int j = 0; j < n_part; j++)
    glvertex3d (view, p.loc[j].x, p.loc[j].y*cos(loc[j].z), p.loc[j].y*sin(loc[j].z));
  glEnd();
  view->ni++; 
}

This is the movie-making event.

event movie( t += 0.05; t <= 100){
  view(phi = 0.5);
  scalar omega[];
  vorticity(u, omega);
  view (fov = 30, quat = {-0.170586,0.0694095,-0.696274,0.693747},
	tx = -0.0209548, ty = -0.204598, width = 650, height = 400);
  draw_vof ("f", lw=6);
  box();
  scatter3Din2D(loc = loc);
  mirror ({-1}){
    draw_vof ("f", lw=6);
    box();
  }
  if (fabs(t-tf) < 1)
    draw_string("Switch BCs", size = 25);
  save ("mov.mp4");
}

Movie

Remember that the flowfield \overrightarrow{u} is 2D: (\overrightarrow{u}(r,z)).

It appears that the tracers are well mixed by the secondary circulation after a few rotations, a nice example of laminar mixing. Consider using this mechanism next time you wish to mix milk in your coffee.