sphere-rotating.c

Sphere rotating in a free-stream flow for Re=200 and \alpha =1

Sphere rotating in a free-stream flow for Re=200 and \alpha =1

This test case is the 3D counterpart of the test case cylinder-rotating.c. We study here the flow induced by the rotation of a sphere in a free-stream flow. Two dimensionless parameters govern this test case:

the Reynolds number Re = u d/\nu;
the non-dimensional rotation rate \alpha = \omega R/u.

We solve here the Navier-Stokes equations and add the sphere using an embedded boundary for Re=200 and \alpha = 1.

#include "grid/octree.h"
#include "../myembed.h"
#include "../mycentered.h"
#include "../myembed-moving.h"
#include "view.h"

Reference solution

#define d    (1.)
#define uref (1.) // Reference velocity, uref
#define tref ((d)/(uref)) // Reference time, tref=d/u
#define Re   (200.) // Particle Reynolds number
#define A    (1.) // Adimensional rotation rate

We also define the shape of the domain.

#define sphere(x,y,z) (sq ((x)) + sq ((y)) + sq ((z)) - sq ((d)/2.))

void p_shape (scalar c, face vector f, coord p)
{
  vertex scalar phi[];
  foreach_vertex()
    phi[] = (sphere ((x - p.x), (y - p.y), (z - p.z)));
  boundary ({phi});
  fractions (phi, c, f);
  fractions_cleanup (c, f,
		     smin = 1.e-14, cmin = 1.e-14);
}

Setup

We need a field for viscosity so that the embedded boundary metric can be taken into account.

face vector muv[];

We define the mesh adaptation parameters.

#define lmin (5) // Min mesh refinement level (l=5 is 2pt/d)
#define lmax (8) // Max mesh refinement level (l=8 is 16pt/d)
#define cmax (1.e-2*(uref)) // Absolute refinement criteria for the velocity field

int main ()
{

The domain is 16\times 16\times 16.

  L0 = 16.;
  size (L0);
  origin (-L0/2., -L0/2., -L0/2.);

We set the maximum timestep.

  DT = 1.e-2*(tref);

We set the tolerance of the Poisson solver.

  TOLERANCE    = 1.e-4;
  TOLERANCE_MU = 1.e-4*(uref);

We initialize the grid.

  N = 1 << (lmin);
  init_grid (N);
  
  run();
}

Boundary conditions

We use inlet boundary conditions.

u.n[left] = dirichlet ((uref));
u.t[left] = dirichlet (0);
u.r[left] = dirichlet (0);
p[left]   = neumann (0);
pf[left]  = neumann (0);

u.n[right] = neumann (0);
u.t[right] = neumann (0);
u.r[right] = neumann (0);
p[right]   = dirichlet (0);
pf[right]  = dirichlet (0);

We give boundary conditions for the face velocity to “potentially” improve the convergence of the multigrid Poisson solver.

uf.n[left]   = (uref);
uf.n[bottom] = 0;
uf.n[top]    = 0;
uf.n[back]   = 0;
uf.n[front]  = 0;

Properties

event properties (i++)
{
  foreach_face()
    muv.x[] = (uref)*(d)/(Re)*fm.x[];
  boundary ((scalar *) {muv});
}

Initial conditions

event init (i = 0)
{

We set the viscosity field in the event properties.

  mu = muv;

We use “third-order” face flux interpolation.

#if ORDER2
  for (scalar s in {u, p, pf})
    s.third = false;
#else
  for (scalar s in {u, p, pf})
    s.third = true;
#endif // ORDER2

We use a slope-limiter to reduce the errors made in small-cells.

#if SLOPELIMITER
  for (scalar s in {u}) {
    s.gradient = minmod2;
  }
#endif // SLOPELIMITER
  
#if TREE

When using TREE and in the presence of embedded boundaries, we should also define the gradient of u at the cell center of cut-cells.

#endif // TREE

We initialize the embedded boundary.

#if TREE

When using TREE, we refine the mesh around the embedded boundary.

  astats ss;
  int ic = 0;
  do {
    ic++;
    p_shape (cs, fs, p_p);
    ss = adapt_wavelet ({cs}, (double[]) {1.e-30},
			maxlevel = (lmax), minlevel = (1));
  } while ((ss.nf || ss.nc) && ic < 100);
#endif // TREE
  
  p_shape (cs, fs, p_p);

We initialize the particle’s rotation speed, in the counter-clockwise direction.

  p_w.x = (A)*(uref)/((d)/2.);
  p_w.z = (A)*(uref)/((d)/2.);

We initialize the velocity.

  foreach()
    u.x[] = cs[]*(uref);    
  boundary ((scalar *) {u});  
}