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. 194). 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. 195 and are in good agreement with the analytical profiles.

Fig. 194 Gravity benchmark for a thick shell. Gravitational potential computed on a line from the center of a constant density shell to a radius of 10,000 km. The gray area indicates the region \(R_{inner} \leq r \leq R_{outer}\) inside the shell, where the density is not zero.

Fig. 195 Gravity benchmark for a thick shell. Gravitational acceleration computed on a line from the center of a constant density shell to a radius of 10,000 km. The gray area indicates the region \(R_{inner} \leq r \leq R_{outer}\) inside the shell, where the density is not zero.