src/examples/yang.c

Flow in a rotating bottom-driven cylindrical container

A cylindrical container is partially filled with liquid. The bottom disk spins at constant speed. A steady state is eventually reached for which the hydrostatic pressure due to the deformation of the free-surface balances the centrifugal force. This is analogous to the solid body rotation in “Newton’s bucket” but more complex due to the development of boundary layers and recirculations generated by the discontinuous velocity boundary condition between the rotating bottom disk and the stationary walls.

This example reproduces one of the cases studied by Wen Yang, 2018 in his PhD.

We use the two-phase axisymmetric Navier–Stokes solver with swirl.

#include "grid/multigrid.h"
#include "axi.h"
#include "navier-stokes/centered.h"
#include "navier-stokes/swirl.h"
#include "two-phase.h"
#include "tension.h"
#include "axistream.h"
#include "view.h"

Boundary conditions

The left boundary is the “bottom” of the container. The azimuthal velocity component w on the bottom disk is \Omega r with the angular velocity \Omega = 1 and r = y.

double Omega = 1.;
u.t[left] = dirichlet(0);
w[left]   = dirichlet(Omega*y);

The “top” boundary is a no-slip stationary wall.

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

The cylinder wall is also no-slip and stationary.

u.t[top] = dirichlet(0);
w[top]   = dirichlet(0);

Parameters

The independent parameters are the aspect ratio of the cylindrical container H/R (unity here), the aspect ratio of the fluid layer \displaystyle G = h/R the Froude number \displaystyle Fr = \frac{R\Omega^2}{g} with g the acceleration of gravity, the Reynolds number of the liquid \displaystyle Re = \frac{R^2\Omega}{\nu_l} the Weber number \displaystyle We = \frac{\rho_l\Omega^2R^3}{\sigma} the ratio of gas to liquid densities and dynamic viscosities \displaystyle \rho_r = \rho_g/\rho_l \displaystyle \mu_r = \mu_g/\mu_l The values correspond to those used by Yang for figure 5.9 in his PhD (see also p. 106).

double R = 1 [1],
       G = 0.25 [0],
       Fr = 0.88 [0],
       Re = 3063 [0],
       We = 3153 [0],
       rhor = 1.205/1.2107e3 [0],
       mur = 18.2e-6/6.09e-2 [0];

Setup

The domain is the square box of unit size (i.e. R = 1), and the angular velocity \Omega is one (see the boundary condition for w above). The gives the values for viscosities, densities and surface tension below. Note that surface tension is commented out since it has very little influence but slows down the calculation.

int main()
{
  N = 128;
  mu1 = sq(R)*Omega/Re;
  rho2 = rho1*rhor;
  mu2 = mu1*mur;
  //  f.sigma = 1./We;
  DT = 0.1 [0,1];
  run();
}

We add gravity (given by the Froude number).

event acceleration (i++)
{
  face vector av = a;
  foreach_face(x)
    av.x[] -= R*sq(Omega)/Fr;
}

The initial interface is flat and positioned at z = RG.

event init (i = 0) {
  fraction (f, R*G - x);
}

#if 0
event snapshots (i += 100) {
  scalar psi[];
  psi[left] = dirichlet (0);
  psi[right] = dirichlet (0);
  psi[top] = dirichlet (0);
  axistream (u, psi);
  dump();
}
#endif

Results

We output the timestep, the total kinetic energy and the total mass.

event logfile (i += 10)
{
  double ke = 0.;
  foreach (reduction(+:ke))
    ke += dv()*rho(f[])*(sq(u.x[]) + sq(u.y[]) + sq(w[]));
  fprintf (stderr, "%g %g %g %g\n", t, dt, ke, statsf (f).sum);
}

We make a movie of the streamlines and isolines of the azimuthal velocity component. They can be directly compared with figure 5.9 of Yang, 2018, reproduced below. Results are close but not identical.

Streamlines and interface reproduced from Yang, 2018, Figure 5.9.

event movie (t += 1; t <= 300)
{
  view (fov = 22.7232, quat = {0,0,-0.707,0.707},
	tx = -0.12824, ty = -0.446981,
	width = 830, height = 432);
  box();

We compute the axisymmetric streamfunction \psi from the velocity components.

  scalar psi[];
  psi[left] = dirichlet (0);
  psi[right] = dirichlet (0);
  psi[top] = dirichlet (0);
  axistream (u, psi);

The values for the isolines are the same as in Yang, 2018, caption of Figure 5.9, p. 107.

  squares ("psi", linear = true, spread = -1);
  isoline ("psi", n = 21, min = 0, max = 0.006);
  isoline ("psi", -0.0012, lc = {0,1,0});
  draw_vof ("f", lc = {1,0,0}, lw = 2);
  translate (y = -1.15) {
    box();
    squares ("w", linear = true, spread = -1);
    isoline ("w", n = 21, min = 0, max = 1);
    draw_vof ("f", lc = {1,0,0}, lw = 2);
  }
  save ("movie.mp4");
}

#if 0 // TREE
event adapt (i++) {
  adapt_wavelet ({u,w}, (double[]){1e-3,1e-3,1e-3}, minlevel = 5, maxlevel = 8);
}
#endif

References