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.
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
[yang2018] |
Wen Yang. Études expérimentales et numériques d’écoulements tournants à surface libre déformable. PhD thesis, Sorbonne Université, December 2018. |