# 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$$:

$\frac{1}{r^2} \frac{\partial }{\partial r} \left(r^2 \frac{\partial U}{\partial r} \right) = 0.$

The solution to this expression is:

(54)#$g=\frac{\partial U}{\partial r} = \frac{C}{r^2},$

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

$g=\frac{\partial U}{\partial r} = \frac{4 \pi}{3} G \rho_0 r + \frac{A}{r^2},$

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$$,

(55)#$g=\frac{\partial U}{\partial r} = \frac{4 \pi}{3} G \rho_0 \left(r - \frac{R_\text{inner}^3}{r^2} \right).$

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}$$:

(56)#$g(r) = \frac{G M}{r^2},$

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)):

$U(r)=-\frac{GM}{r} +D,$

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):

$U(r)= \frac{4 \pi}{3} G \rho_0 \left(\frac{r^2}{2} + \frac{R_\text{inner}^3}{r} \right) + E,$

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

$U(r)=2\pi G \rho_0 (R_\text{inner}^2-R_\text{outer}^2).$

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.