Stokes flow past a periodic array of cylinders

We compare the numerical results with the solution given by the multipole expansion of Sangani and Acrivos, 1982.

#include "embed.h"
#include "navier-stokes/centered.h"
#include "view.h"

This is Table 1 of Sangani and Acrivos, 1982, where the first column is the volume fraction Φ of the cylinders and the second column is the non-dimensional drag force per unit length of cylinder F/(μU) with μ the dynamic vicosity and U the average fluid velocity.

static double sangani[9][2] = {
  {0.05, 15.56},
  {0.10, 24.83},
  {0.20, 51.53},
  {0.30, 102.90},
  {0.40, 217.89},
  {0.50, 532.55},
  {0.60, 1.763e3},
  {0.70, 1.352e4},
  {0.75, 1.263e5}

We will vary the maximum level of refinement, nc is the index of the case in the table above, the radius of the cylinder will be computed using the volume fraction Φ.

int maxlevel = 6, nc;
double radius;

This function defines the embedded volume and face fractions.

void cylinder (scalar cs, face vector fs)
  vertex scalar φ[];
    φ[] = sq(x) + sq(y) - sq(radius);
  fractions (φ, cs, fs);

int main()

The domain is the periodic unit square, centered on the origin.

  size (1.);
  origin (-L0/2., -L0/2.);
  periodic (right);
  periodic (top);

We turn off the advection term. The choice of the maximum timestep and of the tolerance on the Poisson and viscous solves is not trivial. This was adjusted by trial and error to minimize (possibly) splitting errors and optimize convergence speed.

  stokes = true;
  DT = 1e-3;

We do the 9 cases computed by Sangani & Acrivos. The radius is computed from the volume fraction.

  for (nc = 0; nc < 9; nc++) {
    maxlevel = 6;
    N = 1 << maxlevel;
    radius = sqrt(sq(L0)*sangani[nc][0]/π);

We need an extra field to track convergence.

scalar un[];

event init (t = 0)

We initialize the embedded geometry.

  cylinder (cs, fs);

And set acceleration and viscosity to unity.

  const face vector g[] = {1.,0.};
  a = g;
  μ = fm;

The boundary condition is zero velocity on the embedded boundary.

  u.n[embed] = dirichlet_embed(0);
  u.t[embed] = dirichlet_embed(0);

We initialize the reference velocity.

    un[] = u.y[];

We check for a stationary solution.

event logfile (i++; i <= 500)
  double avg = normf(u.x).avg, du = change (u.x, un)/(avg + SEPS);
  fprintf (fout, "%d %d %d %d %d %d %d %d %.3g %.3g %.3g %.3g %.3g\n",
	   maxlevel, i,
	   mgp.i, mgp.nrelax, mgp.minlevel,
	   mgu.i, mgu.nrelax, mgu.minlevel,
	   du, mgp.resa*dt, mgu.resa, statsf(u.x).sum, normf(p).max);

  if (i > 1 && du < 1e-3) {

We output the non-dimensional force per unit length on the cylinder F/(μU), together with the corresponding value from Sangani & Acrivos and the relative error.

    stats s = statsf(u.x);
    double Φ = 1. - s.volume/sq(L0);
    double U = s.sum/s.volume;
    double F = sq(L0)/(1. - Φ);
    fprintf (ferr,
	     "%d %g %g %g %g\n", maxlevel, sangani[nc][0], F/U, sangani[nc][1],
	     fabs(F/U - sangani[nc][1])/sangani[nc][1]);

We dump the simulation and draw the mesh for one of the cases.

    p.nodump = false;
    if (maxlevel == 9 && nc == 8) {
      view (fov = 9.78488, tx = 0.250594, ty = -0.250165);
      draw_vof ("cs", "fs", lc = {1,0,0}, lw = 2);
      squares ("u.x", linear = 1, spread = -1);
      save ("mesh.png");

We stop at level 9.

    if (maxlevel == 9)
      return 1; /* stop */

We refine the converged solution to get the initial guess for the finer level. We also reset the embedded fractions to avoid interpolation errors on the geometry.

#if 0
    refine (level < maxlevel);
    adapt_wavelet ({cs,u}, (double[]){1e-2,2e-6,2e-6}, maxlevel);
    cylinder (cs, fs);

The non-dimensional drag force per unit length closely matches the results of Sangani & Acrivos. For Φ=0.75 and level 8 there is only about 6 grid points in the width of the gap between cylinders.

Non-dimensional drag force per unit length

Non-dimensional drag force per unit length

This can be further quantified by plotting the relative error. It seems that at low volume fractions, the error is independent from the mesh refinement. This may be due to other sources of errors, such as the splitting error in the projection scheme. This needs to be explored further.

Relative error on the drag force

Relative error on the drag force

The adaptive mesh for 9 levels of refinement illustrates the automatic refinement of the narrow gap between cylinders.

Adaptive mesh at level 9 for \Phi=0.75.

Adaptive mesh at level 9 for Φ=0.75.