Stokes flow past a periodic array of spheres

This is the 3D equivalent of flow past a periodic array of cylinders.

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

Note that we do not use an adaptive mesh since the 3D gaps are much wider than for the 2D case.

#include "grid/multigrid3D.h"
#include "embed.h"
#include "navier-stokes/centered.h"

This is adapted from Table 2 of Zick and Homsy, 1982, where the first column is the volume fraction Φ of the spheres and the second column is the drag coefficient K such that the force exerted on each sphere in the array is: F=6πμaKU with a the sphere radius, μ the dynamic vicosity and U the average fluid velocity.

static double zick[7][2] = {
  {0.027,   2.008},
  {0.064,   2.810},
  {0.125,   4.292},
  {0.216,   7.442},
  {0.343,  15.4},
  {0.45,   28.1},
  {0.5236, 42.1}

We can 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 = 5, nc;
double radius;

This function defines the embedded volume and face fractions.

void sphere (scalar cs, face vector fs)
  vertex scalar φ[];
    φ[] = sq(x) + sq(y) + sq(z) - sq(radius);
  boundary ({φ});
  fractions (φ, cs, fs);

int main()

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

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

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 = 2e-2;
  NITERMIN = 10;

We do the 7 cases computed by Zick & Homsy. The radius is computed from the volume fraction.

  for (nc = 0; nc < 7; nc++) {
    maxlevel = 5;
    N = 1 << maxlevel;
    radius = pow (3.*zick[nc][0]/(4.*π), 1./3.);

We need an extra field to track convergence.

scalar un[];

event init (t = 0)

We initialize the embedded geometry.

  sphere (cs, fs);

And set acceleration and viscosity to unity.

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

The boundary condition is zero velocity on the embedded boundary.

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

We initialize the reference velocity.

    un[] = u.x[];

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);
  fflush (fout);
  if (i > 1 && du < 1e-3) {

K is computed using formula 4.2 of Zick an Homsy, although the 1Φ factor is a bit mysterious.

    stats s = statsf(u.x);
    double Φ = 4./3.*π*cube(radius)/cube(L0);
    double V = s.sum/s.volume;
    double dp = 1., μ = 1.;
    double K = dp*sq(L0)/(6.*π*μ*radius*V*(1. - Φ));
    fprintf (ferr,
	     "%d %g %g %g %g %g\n",
	     maxlevel, radius, Φ, V, K,

We stop.

    return 1; /* stop */

The drag coefficient closely matches the results of Zick & Homsy.

Drag coefficient as a function of volume fraction

Drag coefficient as a function of volume fraction

This can be further quantified by plotting the relative error. Better than second-order convergence with spatial resolution is obtained.

Relative error on the drag coefficient

Relative error on the drag coefficient



AS Sangani and A Acrivos. Slow flow through a periodic array of spheres. International Journal of Multiphase Flow, 8(4):343-360, 1982.


A.A. Zick and G.M. Homsy. Stokes flow through periodic arrays of spheres. Journal of Fluid Mechanics, 115(1):13-26, 1982.

See also