# Thick shell gravity benchmark#

*This section was contributed by Cedric Thieulot.*

This benchmark tests the accuracy of the gravity field and gravitational potential computed by the gravity postprocessor inside and outside an Earth-sized planet (without its core) of constant density. The domain is a spherical shell with inner radius \(R_\text{inner}=3840~\text{ km}\) and outer radius \(R_\text{outer}=6371~\text{ km}\). The density is constant in the domain and set to \(\rho_0=3300~\text{ kg}/\text{m}^3\).

First, let us calculate the exact profile which we expect the benchmark to reproduce. The gravitational potential \(U\) of a spherically symmetric object satisfies the Poisson equation \(\Delta U = 4\pi G \rho(\mathbf r)\). For a constant density shell, this equation can be solved analytically for the gravitational acceleration and potential inside and outside the planet. Inside (\(r<R_\text{inner}\)) and outside (\(r>R_\text{outer}\)) the spherical shell (i.e. where \(\rho=0\)) the Poisson equation simplifies to the Laplace equation \(\Delta U=0\):

The solution to this expression is:

where \(C\) is a constant of integration. In order to avoid an infinite gravity field at \(r=0\) (where the density is also zero in this particular setup of a shell), we need to impose \(C=0\), i.e. the gravity is zero for \(r\leq R_\text{inner}\). Another way of arriving at the same conclusion is to realize that \(g\) is zero at the center of the body because the material around it exerts an equal force in every direction. Inside the shell, \(\rho=\rho_0\), yielding

where \(A\) is another integration constant. At the inner boundary, \(r=R_\text{inner}\) and \(g=0\), allowing \(A\) to be computed. Substituting in the value of \(A\),

When \(r\geq R_\text{outer}\), the gravitational potential is given by (54). Requiring the gravity field to be continuous at \(r=R_\text{outer}\):

where \(M=\frac{4 \pi}{3} \rho_0(R_\text{outer}^3-R_\text{inner}^3)\) is the mass contained in the shell. For \(r\ge R_\text{outer}\), the potential is obtained by integrating (56)):

where \(D\) is an integration constant which has to be zero since we require the potential to vanish for \(r\rightarrow \infty\). The potential within the shell, \(R_\text{inner}\leq r \leq R_\text{outer}\), is found by integrating (55):

where \(E\) is a constant. Continuity of the potential at \(r=R_\text{outer}\) requires that \(E=-2\pi\rho_0 G R_\text{outer}^2\). Gravitational acceleration is zero for \(r\leq R_\text{inner}\), so the potential there is constant and a continuity requirement yields

The gravity postprocessor in can be used to calculate the radial components of
gravity (\(g_r\) and \(U\)) at arbitrary points using the sampling scheme
‘*list of points*.’ For this benchmark we calculate points along a
line from the center of the planet to a distant point, \(r=0\) to
\(r=10,000~\text{ km}\) (Fig. 225). Arbitrarily, the latitude and longitude
are both set to \(13\text{ \degree}\) so as to avoid potential measurement
artifacts due to symmetry. The list of radii is defined as follows:

```
subsection Postprocess
set List of postprocessors = gravity calculation,visualization
subsection Gravity calculation
set Sampling scheme = list of points
set List of radius = 0,1e6,2e6,3e6,3.5e6,4e6,4.5e6,5e6,5.5e6,6e6,6.371e6,6.5e6,7e6,8e6,9e6,10e6
set List of longitude = 13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13
set List of latitude = 13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13
set Quadrature degree increase = 1
end
end
```

The resulting measurements obtained for a mesh composed of 12 caps of \(32^3\) cells (i.e., 393,216 total mesh cells) are shown in Fig. 226 and are in good agreement with the analytical profiles.