Here's an asymptotic bound - hopefully it's tight.
We throw $N$ little balls of diameter $D$ randomly (uniformly) inside the big sphere of radius $R$.
UPDATE: The new version (lower half) is better than what follows.
Neglecting border effects (reasonable if $N$ is large) the probability that the ball $i$ is "free" (no other overlaps with it) is
$$P(F_i)=\left(1-v(D)\right)^{N-1} \tag{1}$$
where $v(D) \triangleq D^3/R^3$
The probability that all balls are free can be bounded as
$$P(\cap F_i)=1- P(\cup F_i^c)\ge 1 - N(1-P(F_i)) \triangleq g(D,N) \tag{2} $$
For large $N$
$$g(D,N) \approx 1 - N^2 \, v(D) \tag{3} $$
in the range where this is positive, ie. $0\le D \le D_1 \triangleq R/N^{2/3} $
Now, let $t$ be the minimun distance between the sphere centers. Then
$$P(t \ge D) = P(\cap F_i) \ge g(D,N) \tag{4}$$
And then
$$E(t) = \int_0^{\infty} P(t \ge D) dD \ge\int_0^{D_1} g(D,N) \, dD \approx D_1- N^2 \frac{ D_1^4}{4 R^3} = \frac{3 }{ 4 } \frac{R}{N^{2/3}} $$
(Simulation data suggests that the order is right, and so is the bound, but the real coefficient is around $1.12$ - perhaps $9/8$)
Update: (Improved version)
A better approach can be obtained by considering instead of $F_i$ (free ball) the event $S_j\equiv$ "separated pair" (pair of balls are separated, they do not overlap) where $j$ indexes the $M=N(N-1)/2 \approx N^2/2$ pairs.
By the same reasoning:
$$P(S_j)=1-v(D) =1 - \frac{D^3}{R^3} \tag{5}$$
$$P(\cap S_i) \ge \max(1 - M(1-P(S_i)),0)= \max\left(1 - M \frac{D^3}{R^3},0\right) \triangleq h(D,M) \tag{6} $$
The range where $h(D,M)$ is positive, ie. $0\le D \le D_2 \triangleq R/M^{1/3} $
Now, let $t$ be the minimun distance between the sphere centers. Then
$$P(t \ge D) = P(\cap S_i) \ge h(D,M) \tag{7}$$
And then
$$E(t) = \int_0^{\infty} P(t \ge D) dD \ge\int_0^{D_2} h(D,M) \, dD =\\
= \frac{3}{4} \frac{R}{M^{1/3}}
\approx 0.945 \frac{R}{\sqrt[3]{N(N-1)}} \approx 0.945 \frac{R}{N^{2/3}} \tag{8}$$
Update 2 : A simple heuristic which seems to produce the correct coefficient:
Following the approach above, we could assume that $S_i$ are asympotically independent, and then:
$$P(\cap S_i) \approx \left(1-\frac{D^3}{R^3}\right)^M \tag{9}$$
Then
$$E(t) \approx \int_0^{R}\left(1-\frac{D^3}{R^3}\right)^M dD =\\= R \, \Gamma(4/3) \frac{\Gamma(M+1)}{\Gamma(M+4/3)} \approx R \, \Gamma(4/3) M^{-1/3} \approx 1.12508368 \frac{R}{N^{2/3}} \tag{10}$$
Update 3 : Regarding corrections for border effects.
(Lets assume $R=1$ to save notation, it's just a scale factor)
If we wished to include border effects we should replace $(5)$ (computing the balls intersection as here) by
$$1-D^3+\frac{9}{16}D^4 -\frac{1}{32}D^6 \hspace{1cm} 0\le D\le 2$$
The integral gets more complicated, but the (first order) asymptotic result is not altered:
Lemma: For any positive differentiable function $g(x)$ which , in $[0,+\infty)$, has global maximum at $g(0)=1$, and which has zero first and second derivates $g(x)=1-a x^3 + O(x^4)$ we have (variation of Laplace method, see eg here sec 2.1.3)
$$ \int_0^\infty g(x)^M dx = \frac{\Gamma(1/3)}{3 a^{1/3}} M^{-1/3}+ o(M^{-1/3})$$
which again leads as to $(10)$.
Best Answer
Unfortunately, the answer it not what you would like to hear.
There is nothing more to do to your integral except to express it as the Complete elliptic integral of the second kind.
Using $\cos\phi$ even and substituting $\phi=2t$ you get:
$\int_{0}^{2\pi} \sqrt{r^{2} + R^{2} + 2 r R \cos \phi} \ \mathrm{d} \phi$
$=2\int_{0}^{\pi} \sqrt{r^{2} + R^{2} + 2 r R \cos \phi} \ \mathrm{d} \phi$
$=4\int_{0}^{\pi/2} \sqrt{r^{2} + R^{2} + 2 r R \cos(2t)} \ \mathrm{d} t$
$=4\int_{0}^{\pi/2} \sqrt{r^{2} + R^{2} + 2 r R (1-2\sin^2t)} \ \mathrm{d} t$
$=4\int_{0}^{\pi/2} \sqrt{(r+R)^{2} - 4 r R \sin^2t} \ \mathrm{d} t$
$=4(r+R)\int_{0}^{\pi/2} \sqrt{1 - \frac{4 r R}{(r+R)^{2}} \sin^2t} \ \mathrm{d} t$
$=4(r+R)E\left(\frac{2 \sqrt{r R}}{r+R}\right)$
So your answer is $\frac{2}{\pi}(r+R)E\left(\frac{2 \sqrt{r R}}{r+R}\right)$