Recursive Analytical Formulae of Gravitational Fields and Gradient Tensors for Polyhedral Bodies with Polynomial Density Contrasts of Arbitrary Non-negative Integer Orders

Abstract

Exact computation of the gravitational field and gravitational gradient tensor for a general mass body is a core routine to model the density structure of the Earth. In this study, we report on the existence of closed-form solutions of the gravitational potential, gravitational field and gravitational gradient tensor for a general polyhedral mass body with a polynomial density function of arbitrary non-negative integer orders that can simultaneously vary in both horizontal and vertical directions. Our closed-form solutions of the gravitational potential and the gravitational field are singularity-free, which implies that the observation sites can have arbitrary geometric relationships with polyhedral mass source bodies. However, weak logarithmic singularities exist on the edges of polyhedra for the gravitational gradient tensor. A simple prismatic mass body with polynomial density contrast varying in the vertical direction and a complicated dodecahedral mass body with quartic-order density contrasts were tested to verify the accuracy of the newly derived closed-form solutions. For the gravitational potential, gravitational fields and gradient tensors, our closed-form solutions are in excellent agreement with previously published analytical solutions and Gaussian numerical quadrature solutions.

Appendix: Closed-Form Solution for Volume Integral \(I_v(0,0,0,w)\)

Following our previous work (Chen et al. 2018, equation 24 therein), the closed-form and singularity-free solution for volume integral $I_v(0,0,0,w)\,(w=-1,1,3,\ldots ,\infty )$ can be expressed as:

$$\begin{aligned} {I_v(0,0,0,w)=\iiint \limits _{H}R^{w}{\mathrm{d}}v = -\frac{1}{w+3}\sum _{i=1}^{N}h_i I_s(0,0,0, w),} \end{aligned}$$

(43)

transforming the volume integral $I_v(0,0,0,w)$ into N surface integrals $I_s(0,0,0,w)$, which can be calculated by:

$$\begin{aligned} I_s(0,0,0,w)& {} = \iint \limits _{\partial {H}_{i}}R^{w}{\mathrm{d}}s={\frac{1}{w+2}\sum _{j=1}^{M_i}m_{ij}\int _{C_{ij}} \frac{R^{w+2}}{\tau _{ij}^2}{\mathrm{d}}l-\frac{\beta ({\mathbf {o}}_i)}{w+2}|h_i|^{w+2}} \\& {} = {\frac{1}{w+2}\sum _{j=1}^{M_i}B_{ij}^{w+2} -\frac{\beta ({\mathbf {o}}_i)}{w+2}|h_i|^{w+2},} \end{aligned}$$

(44)

where $\tau _{ij}=|{\mathbf {r}}-{\mathbf {o}}_i|$ is the distance from point ${\mathbf {r}}$ to point ${\mathbf {o}}_i$. $\beta ({\mathbf {o}}_i)$ is the solid angle which can be calculated by the following formula:

$$\begin{aligned} {\beta ({\mathbf {o}}_i) = \sum _{j=1}^{M_i}{\hat{{\mathbf {m}}}}_{ij}\cdot \hat{\varvec{\rho }}_{ij}^{\perp }\left( \arctan {\frac{s_{1ij}}{|m_{ij}|}}-\arctan {\frac{s_{0ij}}{|m_{ij}|}}\right) ,} \end{aligned}$$

(45)

where $\hat{\varvec{\rho }}_{ij}^{\perp }$ is the unit vector from point ${\mathbf {o}}_i$ to point ${\mathbf {r}}_{ij}^{\perp }$ (the projection point of point ${\mathbf {o}}_i$ onto edge $C_{ij}$). Note that $\hat{\varvec{\rho }}_{ij}^{\perp }$ is either parallel or antiparallel to ${{\hat{{\mathbf {m}}}}}_{ij}$. The line integral in Eq. (44) is recursively computed by:

$$\begin{aligned} {B_{ij}^{w+2} = m_{ij}\int _{C_{ij}}\frac{R^{w+2}}{\tau _{ij}^2}{\mathrm{d}}l=m_{ij}L(0,w)+h_i^2 B_{ij}^{w}.} \end{aligned}$$

(46)

The initial line integral of the second term on the right-hand side of Eq. (46) is:

$$\begin{aligned} {B_{ij}^{1} =|h_i|\left( \arctan {\frac{|h_i|s_{1ij}}{m_{ij} R_{1ij}}}-\arctan {\frac{|h_i|s_{0ij}}{m_{ij} R_{0ij}}} \right) +m_{ij}\ln {\frac{s_{1ij}+R_{1ij}}{s_{0ij}+R_{0ij}}},} \end{aligned}$$

(47)

and the closed-form solution for the line integral L(0, w) in the first term on the right-hand side is (Gradshteyn and Ryzhik 2007, equation (2.260)):

$$\begin{aligned} {L(0,w)=\int _{C_{ij}}R^{w}{\mathrm{d}}l = \left. \frac{s_{ij}R^{w}}{w+1}\right| _{s_{0ij}}^{s_{1ij}}+\frac{w}{w+1}\left( h_i^2+m_{ij}^2\right) L(0,w-2),} \end{aligned}$$

(48)

with the initial integral being:

$$\begin{aligned} {L(0,1)=\int _{C_{ij}}R {\mathrm{d}}l = \frac{1}{2}\left[ (h_i^2+m_{ij}^2)\ln {\frac{s_{1ij}+R_{1ij}}{s_{0ij}+R_{0ij}}}+s_{1ij}R_{1ij}-s_{0ij}R_{0ij} \right] .} \end{aligned}$$

(49)

Note that the above closed-form expressions for $I_v(0,0,0,w)$ ($w=-1,1,3,\ldots ,\infty$), which are based on the solutions for $I_s(0,0,0,w)$ and L(0, w), are singularity-free. This means that the scalar components $\phi (p,q,t)$ in Eq. (37) and vector components ${\mathbf {g}}(p,q,t)$ in Eq. (38) can be calculated without singularities.

To calculate the gravitational gradient tensor, the tensor components ${\mathbf {T}}(p,q,t)$ in Eq. (39), including the surface integral $h_i I_s(0,0,0,-3)$, have to be evaluated. For $h_i I_s(0,0,0,-3)$, the closed-form solutions in Eqs. (44) and (46) can be applied, but with the following initial integrals:

$$\begin{aligned} {h_i B_{ij}^{-1} = \cot {\frac{|h_i|s_{1ij}}{m_{ij} R_{1ij}}}-\cot {\frac{|h_i|s_{0ij}}{m_{ij} R_{0ij}}},} \end{aligned}$$

(50)

and

$$\begin{aligned} {L(0,-1)=\int _{C_{ij}}R^{-1} {\mathrm{d}}l = \left\{ {\begin{array}{*{20}{l}} \ln \left( \frac{s_{1ij}+R_{1ij}}{s_{0ij}+R_{0ij}}\right) ,&\quad \text {if}\quad (h_i^2+m_{ij}^2)\ne 0,\\ \left| \ln \frac{s_{1ij}}{s_{0ij}}\right|,& \quad {\text {if}}\quad (h_i^2+m_{ij}^2)=0\quad {\text {and}}\quad {{\mathbf {r}}'} \notin C_{ij}. \end{array}}\right. } \end{aligned}$$

(51)

Weak logarithmic singularities exist in computation of line integral $L(0,-1)$ in Eq. (51) when the observation site locates on the edge of the polyhedral body.