src/test/spheres.c

    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 \Phi of the spheres and the second column is the drag coefficient K such that the force exerted on each sphere in the array is: \displaystyle F = 6\pi\mu a K U with a the sphere radius, \mu 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 \Phi.

    int maxlevel = 5, nc;
    double radius;

    This function defines the embedded volume and face fractions.

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

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

      size (1.);
      origin (-L0/2., -L0/2., -L0/2.);
      foreach_dimension()
        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;
      TOLERANCE = HUGE;
      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.*pi), 1./3.);
        run();
      }
    }

    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;
      mu = 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.

      foreach()
        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 - \Phi factor is a bit mysterious.

        stats s = statsf(u.x);
        double Phi = 4./3.*pi*cube(radius)/cube(L0);
        double V = s.sum/s.volume;
        double dp = 1., mu = 1.;
        double K = dp*sq(L0)/(6.*pi*mu*radius*V*(1. - Phi));
        fprintf (ferr,
    	     "%d %g %g %g %g %g\n",
    	     maxlevel, radius, Phi, V, K,
    	     zick[nc][1]);

    We stop.

        return 1; /* stop */
      }
    }

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

    set xlabel 'Volume fraction'
    set ylabel 'K'
    set logscale y
    set grid
    set key top left
    plot 'log' u 3:6 ps 1 lw 2 t 'Zick and Homsy, 1982',	     \
         '' u 3:5 ps 1 pt 6 lw 2 t '5 levels',		     \
         'spheres.6' u 3:5 ps 1 pt 8 lw 2 t '6 levels'
    Drag coefficient as a function of volume fraction (script)

    Drag coefficient as a function of volume fraction (script)

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

    set ylabel 'Relative error'
    plot 'log' u 3:(abs($6-$5)/$5) w lp t '5 levels',	\
         'spheres.6' u 3:(abs($6-$5)/$5) w lp t '6 levels'
    Relative error on the drag coefficient (script)

    Relative error on the drag coefficient (script)

    References

    [sangani1982]

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

    [zick1982]

    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