Falling Sphere Benchmark

This benchmark simulates a rigid sphere falling through a cylinder filled with a viscous medium as in Figure 1.

Figure 1: Schematic of a Sphere falling through a Cylinder.

The file input/benchmarks/falling_sphere/README has instructions on running this benchmark. In an infinitely large cylinder, the analytic solution for the drag on a sphere is

[;F=6\pi\eta r u,;]

where [;\eta;] is the viscosity of the medium, [;r;] is the radius of the sphere, and [;u;] is the velocity of the sphere. Conversely, the buoyancy force is

[;F=\frac{4}{3}\pi{r^3}g\delta\rho,;]

where [;g;] is the gravitational constant and [;\delta\rho;] is the density difference between the sphere and the medium. Balancing these two forces and solving for the velocity gives

[;u = \frac{2}{9}{r^2}g\delta\rho / \eta.;]

Setting [;g=1;], [;r=1;], [;\delta\rho=0.01;], and [;\eta=1;] gives a velocity of

[;u=0.00222.;]

In our case, we simulate a rigid sphere with a high viscosity sphere. This allows some internal circulation within the sphere, and so the expression for the velocity becomes,

[;u = \frac{1}{3}\frac{r^2{g}\delta\rho}{\eta}\frac{\eta+\eta'}{\eta+\frac{3}{2}\eta'},;]

where [;\eta';] is the viscosity of the sphere. For our case, the background medium's viscosity is 1 and the sphere's viscosity is 100, so the correction is about 1%. This turns out to be smaller than other effects for the cases we ran.

When the boundaries are not infinitely far away, we can expand the solution in terms of the ratio of the radius of the sphere ([;r;]) to the radius of the cylinder ([;R;]). One solution by Habermann gives a drag force of

[;F_H=6\pi\eta{ru}\frac{1-0.75857\cdot\left(\frac{r}{R}\right)^5}{1+f_H\left(\frac{r}{R}\right)},;]

where

[;f_H\left(\frac{r}{R}\right) = -2.1050(r/R) + 2.0865(r/R)^3 -1.7068(r/R)^5+0.72603(r/R)^6.;]

For our case with [;r=1;], [;R=4;], this gives a velocity of

[;u=1.122747319\cdot{10^{-3}}.;]

The walls reduce the speed by about a factor of two.

Another solution by Faxen gives a drag force of

[;F_F=6\pi\eta{ru}/(1+f_F(r/R)),;]

where

[;f_f(r/R)=-2.10444(r/R)+2.08877(r/R)^3-0.94813(r/R)^5 \&-&1.372(r/R)^6+3.87(r/R)^8-4.19(r/R)^{10}.;]

For our case, this gives a speed of

[;u=1.12293603939\cdot{10^{-3}},;]

which agrees closely with the result from Habermann.

Another possible artifact is that we do not simulate an infinitely long cylinder. This turns out to be a small effect. We use a cylinder with a height of 8, and place the sphere halfway down. We did runs where the cylinder was twice as tall, and the results were essentially unchanged.

The errors in the computed velocity compared to the Faxen solution are plotted in Figure 2. These were done with resolutions of 8 × 16 × 8, 16 × 32 × 16, 32 × 64 × 32, and 64 × 128 × 64, corresponding to grid sizes ([;h;]) of 0.5, 0.25, 0.125, and 0.0625. Because of the symmetries of the problem we only have to simulate a quarter of the domain. As we increase the resolution (decrease [;h;]), the error decreases. Since we are simulating a high viscosity sphere rather than a completely rigid sphere, the velocity inside the sphere is not uniform. The error bars indicate the variation in velocity across the sphere.

Figure 2: Error in computed velocity vs. resolution.

Scaling the error with resolution gives Figure 3. The error scales linearly with resolution, giving us confidence that we can accurately solve this problem.

Figure 3: As in Figure 2, but with the error scaled with [;h;]. So the higher resolution errors are multiplied by 2, 4 and 8.