Mathematica integral wrong

Not an answer, but this won't fit in comment:

    Integrate[1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2], 
 Element[{x, y, z}, 
  RegionDifference[Ball[{0, 0, 0}, bigR], Ball[{0, 0, 0}, littleR]]], 
 Element[{x1, y1, z1}, Ball[{0, 0, 0}, littleR]], 
 Assumptions -> 0 < littleR < bigR]  (*gives 0*)

and

Integrate[1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2], 
 Element[{x, y, z}, 
  RegionDifference[Ball[{0, 0, 0}, bigR], 
   Ball[{0, 0, 0}, littleR - epsilon]]], 
 Element[{x1, y1, z1}, Ball[{0, 0, 0}, littleR]], 
 Assumptions -> 0 < littleR < bigR && 0 < epsilon < littleR] (*gives 0*)

and

NIntegrate[1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2], 
 Element[{x, y, z}, 
  RegionDifference[Ball[{0, 0, 0}, 2], Ball[{0, 0, 0}, 1 - 0.1]]],
 Element[{x1, y1, z1}, Ball[{0, 0, 0}, 1 - 0.1]]
 ] (*gives 61.2 but with warnings about accuracy*)

I think that 0 is an incorrect result and worth a bug report.

NIntegrate[1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2], 
 Element[{x, y, z}, Ball[{0, 0, 0}, 1]], 
 Element[{x1, y1, z1}, Ball[{0, 0, 0}, 2]]] (*gives 71.6 with warnings*)

but I don't trust that result at all.

The result is totally wrong. This integral isn't too hard to do in spherical coordinates.

Start with the interaction between two spherical shells of radii $r$ and $r_1$. Let $\{x,y,z\}=\{0,0,r\}$ (we pick the $z$-axis such that it contains the point on the first spherical shell) and $\{x_1,y_1,z_1\}=r_1\{\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta\}$. The integrand becomes

1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2] /.
   {x -> 0, y -> 0, z -> r,
    x1 -> r1 Sin[θ] Cos[φ], y1 -> r1 Sin[θ] Sin[φ], z1 -> r1 Cos[θ]} // FullSimplify
(*    1/Sqrt[r^2 + r1^2 - 2 r r1 Cos[θ]]    *)

We have a Jacobian of $4\pi r^2$ for the point on the first sphere (the point on the first sphere can be chosen from anywhere on the ball's surface) and $2\pi r_1^2\sin\theta$ for the point on the second sphere. The integral is therefore, depending on the ordering of the radii,

Assuming[0 < r < r1,
   Integrate[((4 π r^2) (2 π r1^2 Sin[θ]))/Sqrt[r^2 + r1^2 - 2 r r1 Cos[θ]],
             {θ, 0, π}]]
(*    16 π^2 r^2 r1    *)

Assuming[0 < r1 < r,
   Integrate[((4 π r^2) (2 π r1^2 Sin[θ]))/Sqrt[r^2 + r1^2 - 2 r r1 Cos[θ]],
             {θ, 0, π}]]
(*    16 π^2 r r1^2    *)

Together, this result is the shell–shell interaction:

shell[r_, r1_] = 16 π^2 r r1 Min[r, r1]

You see that this interaction is not singular when $r=r_1$, which helps with the integration over the balls:

Assuming[0 < R < R1,
   Integrate[shell[r, r1], {r, 0, R}, {r1, 0, R1}]] // FullSimplify
(*    -(8/15) π^2 R^3 (R^2 - 5 R1^2)    *)

Assuming[0 < R1 < R,
   Integrate[shell[r, r1], {r, 0, R}, {r1, 0, R1}]] // FullSimplify
(*    -(8/15) π^2 R1^3 (-5 R^2 + R1^2)    *)

Together, the ball–ball interaction is therefore

ball[R_, R1_] = 8/15 π^2 Min[R, R1]^3 (5 (R^2 + R1^2) - 6 Min[R, R1]^2);

So, definitely not zero.

Let's try it out numerically:

With[{R = 0.3, R1 = 0.7},
  {NIntegrate[1/Sqrt[(x - x1)^2 + (y - y1)^2 + (z - z1)^2],
              Element[{x, y, z}, Ball[{0, 0, 0}, R]], 
              Element[{x1, y1, z1}, Ball[{0, 0, 0}, R1]]],
   ball[R, R1]}]

(* ...lots of complaining about accuracy... *)

(*    {0.336256, 0.335409}    *)

seems to work fine.

The result is mathematically correct, therefore there is no bug in Mathematica.

What you are trying to compute is this integral $$ \mathcal{I}=\int_{B_R^3}\int_{B_{R_1}^3}\frac{1}{|\mathbf{x}-\mathbf{y}|}d\mathbf{x}^3d\mathbf{y}^3 $$ We can use Legendre expansion that reads $$ \frac{1}{|\mathbf{x}-\mathbf{y}|}=\sum\limits_{n=0}^{\infty}\frac{r_{<}^{n}}{r_{>}^{n+1}}P_n(\cos(\gamma)) $$ where $r_{<}$ is the smaller of $|\mathbf{x}|$ and $|\mathbf{y}|$ and $r_{>}$ is the larger.

Let us choose a fixed spherical coordinate system for $\mathbf{x}$ such that $\theta,\phi$ are the relevant spherical angles (and $r_x$ is the norm). In comparison, we will choose a relative spherical coordinate system for $\mathbf{y}$ such that its $z-$axis will always point towards $\mathbf{x}$: the relevant angles will then be $\gamma$ (angle between these vectors) and $\rho$ (and $r_y$ is the norm). In this coordinate system, we then have $$ \mathcal{I}=\sum\limits_{n=0}^{\infty}\int\limits_{0}^{R}r_x^2dr_x\int\limits_{0}^{R_1}r_y^2dr_y\int\limits_{0}^\pi\sin(\theta)d\theta\int\limits_{0}^{2\pi}d\phi\int\limits_{0}^\pi\sin(\gamma)d\gamma\int\limits_{0}^{2\pi}d\rho\frac{r_{<}^{n}}{r_{>}^{n+1}}P_n(\cos(\gamma)) $$ which is basically $$ \mathcal{I}=\sum\limits_{n=0}^{\infty}\left(\text{gamma independent terms}\right)\int\limits_{0}^\pi\sin(\gamma)P_n(\cos(\gamma))d\gamma $$

This integral is zero. You can also check it explicitly for various $n$, or directly using Mathematica as well:

FullSimplify[Integrate[Sin[\[Gamma]] LegendreP[n, Cos[\[Gamma]]], {\[Gamma], 0, Pi}, Assumptions -> Element[n, NonNegativeIntegers]],  Assumptions -> Element[n, NonNegativeIntegers]]

Therefore, the given integral $\mathcal{I}$ is zero.