Integral $\int_{d_1}^{d_2} \int_{-L/2}^{L/2} \int_{-L/2}^{L/2} \frac{1}{(x^2+y^2+z^2)^3} dx dy dz$

Given some conditions, the analytic form is found:

Assuming[L > 0 && 0 < d1 < d2, 
  Integrate[1/(x^2 + y^2 + z^2)^3,
            {x, d1, d2}, {y, -L/2, L/2}, {z, -L/2, L/2}] // FullSimplify]

$$ \frac{2 \left(\frac{\left(16 d_1^4+2 d_1^2 L^2+L^4\right) \tan ^{-1}\left(\frac{L}{\sqrt{4 d_1^2+L^2}}\right)}{d_1^3 \sqrt{4 d_1^2+L^2}}+L \left(\frac{1}{d_1}-\frac{1}{d_2}\right)+5 \sqrt{2} \tan ^{-1}\left(\frac{\sqrt{2} d_1}{L}\right)-\frac{\left(16 d_2^4+2 d_2^2 L^2+L^4\right) \tan ^{-1}\left(\frac{L}{\sqrt{4 d_2^2+L^2}}\right)}{d_2^3 \sqrt{4 d_2^2+L^2}}-5 \sqrt{2} \tan ^{-1}\left(\frac{\sqrt{2} d_2}{L}\right)\right)}{3 L^3} $$

A similar formula is achieved for $d_1<d_2<0$. As Henrik comments, the integral diverges if $d_1\le0\le d_2$.

Mathematica is able to solve the indefinite integral:

ClearAll[expr];
expr[x_,y_,z_]=Integrate[1/(x^2+y^2+z^2)^3,x,y,z];

Using the indefinite integral expr found above to find the definite integral as follows:

expr[y_,z_]=expr[L/2,y,z]-expr[-L/2,y,z];
expr[z_]=expr[L/2,z]-expr[-L/2,z];
expr[Subscript[d, 2]]-expr[Subscript[d, 1]]//FullSimplify

$\frac{2 \left(\frac{L^2 \sqrt{4 d_1^2+L^2} \tan ^{-1}\left(\frac{L}{\sqrt{4 d_1^2+L^2}}\right)}{d_1^3}+\frac{L-2 \sqrt{4 d_1^2+L^2} \tan ^{-1}\left(\frac{L}{\sqrt{4 d_1^2+L^2}}\right)}{d_1}+\frac{24 d_1 \tan ^{-1}\left(\frac{L}{\sqrt{4 d_1^2+L^2}}\right)}{\sqrt{4 d_1^2+L^2}}+\frac{L^4 \left(-\tan ^{-1}\left(\frac{L}{\sqrt{4 d_2^2+L^2}}\right)\right)-d_2^2 \left(L \sqrt{4 d_2^2+L^2}+2 \left(8 d_2^2+L^2\right) \tan ^{-1}\left(\frac{L}{\sqrt{4 d_2^2+L^2}}\right)\right)}{d_2^3 \sqrt{4 d_2^2+L^2}}+5 \sqrt{2} \left(\tan ^{-1}\left(\frac{\sqrt{2} d_1}{L}\right)-\tan ^{-1}\left(\frac{\sqrt{2} d_2}{L}\right)\right)\right)}{3 L^3}$

---

Now to numerically estimate it, it must be noted that there is a very fast singularity in the middle, i.e. (x,y,z)=(0,0,0).

Plot[1/(x^2)^3,{x,-10,10}]
Plot3D[1/(x^2+y^2)^3,{x,-10,10},{y,-10,10}]

It just continues the same way in higher dimension. So as long as you avoid the point in the middle by some margin it will be able to numerically estimate its value.

NIntegrate[
    1/(x^2+y^2+z^2)^3,
    {x,-1,1},
    {y,-1,1},
    {z,0,2} (*failure expected since it touches the singularity*)
]

fails to converge and very big values overflow the machine number.

NIntegrate[
    1/(x^2+y^2+z^2)^3,
    {x,-1,1},
    {y,-1,1},
    {z,0.1,2} (*big value expected since close to singularity*)
]

522.763

NIntegrate[
    1/(x^2+y^2+z^2)^3,
    {x,-1,1},
    {y,-1,1},
    {z,1,2} (*small value expected since it has very thin tails away from origin / singularity*)
]

0.310657

INTRODUCTION

In order to compute our integral we have to provide appropriate assumptions in adequate functions and so this turns us to another issue stated by Szabolcs in his interesting discussion When and why are Assuming and Assumptions not equivalent? clarified elswhere by Daniel Lichtblau. Therefore our integral provides an important example to the mentioned issue related to advantages of Assumptions in Integrate over the construction Assuming[ ..., Integrate[...]] and vice versa. And so a straightforward and natural approach follows along this way:

Integrate[1/(x^2 + y^2 + z^2)^3, {x, -(L/2), L/2}, {y, -(L/2), L/  2},
                                 {z, d1, d2}, Assumptions -> L > 0 && 0 < d1 < d2]

nontheless I've had to stop this computation since it hasn't been completed after 30 minutes. Therefore our modiffied approach involves a computation in a few steps.

SOLUTION

int1[L_,z_]= Integrate[ 1/(x^2 + y^2 + z^2)^3, {x, -(L/2), L/2}, {y, -(L/2), L/2}, 
                        Assumptions -> L > 0 && z > 0];
int[L_,d1_,d2_] = Integrate[ int1[L, z], {z, d1, d2}, 
                             Assumptions -> 0 < d1 < d2 && L > 0];
integral[L_, d1_, d2_] = 
  FullSimplify[int[L, d1, d2], Assumptions -> L > 0 && 0 < d1 < d2];

i.e.

TraditionalForm[ integral[L, d1, d2]]

It took on my laptop (i3 CPU 1.9 GHz, 4 GB RAM, Windows 10 x 64,Mathematica 11.2)

AbsoluteTiming[ int1[L,z]; int[L,d1,d2]; integral[L,d1,d2]]

{ 316.34, ...}

roughly the same as Roman's approach, which took { 322.287, ...}, howerver my solution provides an alternative approach which avoids possible issues with generic results as mentioned in introduction.