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 Ωr with the angular velocity Ω=1 and r=y.

u.t[left] = dirichlet(0);
w[left]   = dirichlet(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 G=h/R the Froude number Fr=RΩ2g with g the acceleration of gravity, the Reynolds number of the liquid Re=R2Ωνl the Weber number We=ρlΩ2R3σ the ratio of gas to liquid densities and dynamic viscosities ρr=ρg/ρl μr=μg/μl The values correspond to those used by Yang for figure 5.9 in his PhD (see also p. 106).

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

Setup

The domain is the square box of unit size (i.e. R=1), and the angular velocity Ω 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 = 1./Re;
  rho2 = rho1*rhor;
  mu2 = mu1*mur;
  //  f.sigma = 1./We;
  DT = 0.1;
  run();
}

We add gravity (given by the Froude number).

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

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

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

#if 0
event snapshots (i += 100) {
  scalar ψ[];
  ψ[left] = dirichlet (0);
  ψ[right] = dirichlet (0);
  ψ[top] = dirichlet (0);
  axistream (u, ψ);
  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()*ρ(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 (left) and azimuthal velocity (right). Interface in red.

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

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 ψ from the velocity components.

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

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

[yang2018]

Wen Yang. Études expérimentales et numériques d’écoulements tournants à surface libre déformable. PhD thesis, Sorbonne Université, December 2018.