# The "Stokes' law" benchmark
*This section was contributed by Juliane Dannberg.*
Stokes' law was derived by George Gabriel Stokes in 1851 and describes
the frictional force a sphere with a density different than the surrounding
fluid experiences in a laminar flowing viscous medium. A setup for testing
this law is a sphere with the radius $r$ falling in a highly viscous fluid
with lower density. Due to its higher density the sphere is accelerated by the
gravitational force. While the frictional force increases with the velocity of
the falling particle, the buoyancy force remains constant. Thus, after some
time the forces will be balanced and the settling velocity of the sphere $v_s$
will remain constant:
```{math}
:label: eq:stokes-law
\begin{aligned}
\underbrace{6 \pi \, \eta \, r \, v_s}_{\text{frictional force}} =
\underbrace{4/3 \pi \, r^3 \, \Delta\rho \, g,}_{\text{buoyancy force}}
\end{aligned}
```
where $\eta$ is the dynamic viscosity of the fluid, $\Delta\rho$ is the
density difference between sphere and fluid and $g$ the gravitational
acceleration. The resulting settling velocity is then given by
```{math}
:label: eq:stokes-velo
\begin{aligned}
v_s = \frac{2}{9} \frac{\Delta\rho \, r^2 \, g}{\eta}.
\end{aligned}
```
Because we do not take into account inertia in our numerical computation, the
falling particle will reach the constant settling velocity right after the
first timestep.
For the setup of this benchmark, we chose the following parameters:
```{math}
:label: eq:stokes-parameters
\begin{aligned}
r &= 200 \, \text{ km}\\
\Delta\rho &= 100 \, \text{ kg/m}^3\\
\eta &= 10^{22} \, \text{ Pa.s}\\
g &= 9.81 \, \text{ m/s}^2.
\end{aligned}
```
With these values, the exact value of sinking velocity is $v_s =
8.72\times 10^{-10} \, \text{ m/s}$.
To run this benchmark, we need to set up an input file that describes the
situation. In principle, what we need to do is to describe a spherical object
with a density that is larger than the surrounding material. There are
multiple ways of doing this. For example, we could simply set the initial
temperature of the material in the sphere to a lower value, yielding a higher
density with any of the common material models. Or, we could use ASPECT's
facilities to advect along what are called "compositional fields"
and make the density dependent on these fields.
We will go with the second approach and tell to advect a single compositional
field. The initial conditions for this field will be zero outside the sphere
and one inside. We then need to also tell the material model to increase the
density by $\Delta\rho=100 \text{ kg m}^{-3}$ times the concentration of the
compositional field. This can be done, like everything else, from the input
file.
All of this setup is then described by the following input file. (You can find
the input file to run this cookbook example in
[cookbooks/stokes/stokes.prm](https://www.github.com/geodynamics/aspect/blob/main/cookbooks/stokes/stokes.prm).
For your first runs you will probably want to
reduce the number of mesh refinement steps to make things run more quickly.)
```{literalinclude} stokeslaw.prm
```
Using this input file, let us try to evaluate the results of the current
computations for the settling velocity of the sphere. You can visualize the
output in different ways, one of it being ParaView and shown in
{numref}`fig:stokes-falling-sphere-2d` (an alternative is to use Visit as described in
{ref}`sec:run-aspect:visualizing-results`; 3d images of this simulation using Visit are
shown in {numref}`fig:stokes-falling-sphere-3d`). Here, ParaView has the advantage that you can
calculate the average velocity of the sphere using the following filters:
1. Threshold (Scalars: C_1, Lower Threshold 0.5, Upper Threshold 1),
2. Integrate Variables,
3. Cell Data to Point Data,
4. Calculator (use the formula sqrt(velocity_x^2+
velocity_y^2+velocity_z^2)/Volume).
If you then look at the Calculator object in the Spreadsheet View, you can see
the average sinking velocity of the sphere in the column "Result"
and compare it to the theoretical value
$v_s = 8.72\times 10^{-10} \, \text{ m/s}$. In this case, the numerical result is
$8.865\times 10^{-10} \,
\text{ m/s}$ when you add a few more refinement steps to actually resolve
the 3d flow field adequately. The "velocity statistics"
postprocessor we have selected above also provides us with a maximal velocity
that is on the same order of magnitude. The difference between the analytical
and the numerical values can be explained by different at least the following
three points: (i) In our case the sphere is viscous and not rigid as assumed
in Stokes' initial model, leading to a velocity field that varies inside
the sphere rather than being constant. (ii) Stokes' law is derived using
an infinite domain but we have a finite box instead. (iii) The mesh may not
yet fine enough to provide a fully converges solution. Nevertheless, the fact
that we get a result that is accurate to less than 2% is a good indication
that implements the equations correctly.
```{figure-md} fig:stokes-falling-sphere-2d
Stokes benchmark. Both figures show only a 2D slice of the three-dimensional model. Left: The compositional field and overlaid to it some velocity vectors. The composition is 1 inside a sphere with the radius of 200 km and 0 outside of this sphere. As the velocity vectors show, the sphere sinks in the viscous medium. Right: The density distribution of the model. The compositional density contrast of 100 kg/\text{ m}^3 leads to a higher density inside of the sphere.
```
```{figure-md} fig:stokes-falling-sphere-3d
Stokes benchmark. Three-dimensional views of the compositional field (left), the adaptively refined mesh (center) and the resulting velocity field (right).
```