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
As Parcly Taxel pointed out, MathWorld has a page on Hypercube Line Picking, with many references.
Mathworld gives a table of the mean distance for hypercubes up to $n=8$ dimensions. For reasons I will explain below, for large $n$ a good approximation of the mean is $\sqrt{\dfrac{n}{6}-\dfrac{7}{120}}$ and it is not bad for small $n$ either
For $n=1$, you have a triangular distribution for the distance with density $f_{d_1}(x)=2-2d_1$ for $0 \lt x \le 1$, giving a mean of $\frac13$, a variance of $\frac1{18}$ and a second moment of $\frac1{6}$. The square of the distance has density $f_{d_1^2}(x)=\frac{1}{\sqrt{x}}-1$ for $0 \lt x\le 1$, giving a mean of $\frac16$, a variance of $\frac7{180}$ and a second moment of $\frac1{15}$.
It gets more complicated for higher dimensions, but (as Ivan Neretin says) the Central Limit Theorem tells us that the square of the distance is almost normally distributed for large $n$, with mean $\frac{n}{6}$ and variance $\frac{7n}{180}$. So we can say $$\dfrac{D_n^2 - \frac{n}{6}}{\sqrt{\frac{7n}{180}}} \ \xrightarrow{d}\ N(0,1)$$
Less obviously, the distance itself is also almost normally distributed for large $n$. In general we can say that if $X_1, \ldots, X_n$ are i.i.d. random variables with finite non-zero mean $\mu$ and variance $\sigma^2$, and $\displaystyle Y=\sum_{i=1}^n X_i$ and $Z=\sqrt{|Y|}$, then $\displaystyle \dfrac{Z - \sqrt{n |\mu|-\tfrac{\sigma^2}{4|\mu|}}}{\sqrt{\tfrac{\sigma^2}{4|\mu|}}}\ \xrightarrow{d}\ N(0,1)$ as $n$ increases. In this particular case $\mu=\frac{1}{6}$ and $\sigma^2 = \frac7{180}$ as statistics of the $1$-dimensional square of distance, so we can say $$\dfrac{D_n - \sqrt{\frac{n}{6}-\frac7{120}}}{\sqrt{\frac{7}{120}}} \ \xrightarrow{d}\ N(0,1)$$ suggesting an approximate mean for the distance of $\sqrt{\frac{n}{6}-\frac7{120}}$ and approximate variance of $\frac7{120}$ when $n$ is large.
The actual densities are not simple analytically, but the following graph uses numerical convolution and integration to illustrate the densities for the distance when $n=1$ to $16$ and also shows in red the normal approximation when $n=16$.
For large $n$ the variance of the distance stays close to $\frac{7}{120}$ making the standard deviation about $0.24$. For example, with an $n=2500$ dimensional unit hypercube, the distance can be anything from $0$ to $50$ but in the large majority of cases it will be between $20$ and $21$ and in all but a vanishingly tiny proportion of cases it will be between $19$ and $22$. In data analysis, this curse of dimensionality means there can be relatively little difference in the distances between different pairs of random samples.