sphere-oscillating.c

Flow past a oscillating sphere for Re=40 and St=3.2

Flow past a oscillating sphere for Re=40 and St=3.2

This test case is based on the numerical work of Mei, 1994 and is the 3D counterpart of the test case cylinder-oscillating.c. We reproduce here the results of Mei, 1994, Wang et al., 2010 and Schneiders et al., 2013. The authors studied the flow induced by the sinusoidal oscillation of a sphere in the x-direction.

We solve here the Navier-Stokes equations and add the sphere using an embedded boundary.

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

Reference solution

#define d    (1.)
#define Re   (40.) // Particle Reynolds number
#define p_w0 (3.2) // Particle pulsation
#define uref (1.) // Reference velocity, uref
#define tref (min ((d)/(uref), 2.*M_PI/(p_w0))) // Reference time, tref=min (d/u, 2pi/w)

We define the sphere’s imposed motion.

#define p_acceleration(t,w) ((w)*cos((w)*(t))) // Particle acceleration
#define p_velocity(t,w)     (sin((w)*(t))) // Particle velocity
#define p_displacement(t,w) (-1./(w)*cos((w)*(t))) // Particle displacement

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 and vary the maximum level of refinement.

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

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-6;
  TOLERANCE_MU = 1.e-6*(uref);  
  
  for (lmax = 8; lmax <= 10; lmax++) {

We initialize the grid.

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

Boundary conditions

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

uf.n[left]   = 0;
uf.n[right]  = 0;
uf.n[bottom] = 0;
uf.n[top]    = 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.

We first define the sphere’s initial position.

  p_p.x = (p_displacement(0., (p_w0)));
  
#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 speed and acceleration.

  p_au.x = (p_acceleration (0., (p_w0)));
  p_u.x  = (p_velocity     (0., (p_w0)));
}

Embedded boundaries

The particle’s position is advanced to time t + \Delta t.

event advection_term (i++)
{
  p_au.x = (p_acceleration (t + dt, (p_w0)));
  p_u.x  = (p_velocity     (t + dt, (p_w0)));
  p_p.x  = (p_displacement (t + dt, (p_w0)));
}

We verify here that the velocity and pressure gradient boundary conditions are correctly computed.

event check (i++)
{
  foreach() {
    if (cs[] > 0. && cs[] < 1.) {

      // Normal pointing from fluid to solid
      coord b, n;
      embed_geometry (point, &b, &n);

      // Velocity
      bool dirichlet;
      double ub;