# 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)).

The bounding domain (box), the free surface (black line) and flow tracers (dots)

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.