This is very classical (also known as "Thomson problem" usually on the sphere but the same asymptotic results hold for other sets as well - Hausdorff dimension is important).

Whenever

$$
\inf_{\mu}\int \int_{A\times A}\frac{1}{|x-y|^{s}} d\mu(x) d\mu(y)<\infty,
$$

where infimum is taken with respect to all probability measures supported on, say a compact set $A \subset \mathbb{R}^{k}$, (this usually happens when $0\leq s <d_{H}(A),$ here $d_{H}(A)$ is the Hausdorff dimension of the set $A$) then one has
$\frac{1}{N^{2}} \inf_{x_{1}, \ldots, x_{N} \in A} \sum_{1\leq i \neq j \leq N}\frac{1}{|x_{i}-x_{j}|^{s}} \to \inf_{\mu}\int_{A\times A}\frac{1}{|x-y|^{s}} d\mu(x) d\mu(y)$
as $N$ goes to infinity. In other words $s$-Riesz energy
$$
E_{s}(A):=\inf_{x_{1}, \ldots, x_{N} \in A} \sum_{1\leq i \neq j \leq N}\frac{1}{|x_{i}-x_{j}|^{s}} \sim C N^{2}.
$$

The standard references is potential theory

*Landkof, N. S.*, Foundations of modern potential theory. Translated from the Russian by A. P. Doohovskoy, Die Grundlehren der mathematischen Wissenschaften. Band 180. Berlin-Heidelberg-New York: Springer-Verlag. X, 424 p. Cloth DM 88.00; $ 27.90 (1972). ZBL0253.31001., Springer-Verlag

*Mattila, Pertti*, Geometry of sets and measures in Euclidean spaces. Fractals and rectifiability, Cambridge Studies in Advanced Mathematics. 44. Cambridge: Univ. Press. xii, 343 p. (1995). ZBL0819.28004.

In your case $s=1$ and $d_{H}([0,1]^{d})=d$.

Another interesting scenario is when $s=d$ (in this case $E_{s}(A)\sim C n^{2} \ln n$); and when $s>d$ (in this case $E_{s}(A) \sim C N^{1+\frac{s}{d}}$). You may look at

*Hardin, D. P.; Saff, E. B.*, **Discretizing manifolds via minimum energy points**, Notices Am. Math. Soc. 51, No. 10, 1186-1194 (2004). ZBL1095.49031.

and references therein.

1more comment