/** ![Some people prefer to mix some milk in their coffee. Image via [Antony Hall](http://www.antonyhall.net/).](https://i1.wp.com/antonyhall.net/images/coffee-cup-2015/4.jpg) # 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: $$\frac{h_0}{R}= 0.5$$ $$\frac{\rho_1}{\rho_2} = 1000$$ $$\frac{\mu_1}{\mu_2} = 50$$ $$\frac{\omega^2 R}{g} = 0.2$$ $$\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)](coffee/mov.mp4) 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. */