Minimal distance limit problem
The limit does exist and is equal to $1/4$. We can prove this with a bit of equidistribution theory. Given a doubly indexed sequence $x_{j,n}$ of real numbers over $1\le j \le n$, we shall say that this is equidistributed (mod 1) as $n\to\infty$ if $$ \lim_{n\to\infty}\frac1n\sum_{j=1}^nF(x_{j,n})=\int_0^1F(x)\,dx $$ for all continuous $F\colon\mathbb{R}\to\mathbb{C}$ with period 1. In particular, the function $F(x)=\min_{k\in\mathbb{Z}}\lvert x-k\rvert$ satisfies the required properties so, with $d_j$ as in the question, $$ \sum_{j=1}^nd_j = \sum_{j=1}^n\frac1nF\left(n\sqrt{1-(j/n)^2}\right)=\frac1n\sum_{j=1}^nF\left(\sqrt{n^2-j^2}\right). $$ Therefore, this converges to $\int_0^1F(x)\,dx=1/4$ as long as the sequence $x_{j,n}=\sqrt{n^2-j^2}$ ($1\le j\le n$) is equidistributed mod 1 as $n\to\infty$. I'll show this using the following statement.
Theorem: Let $f\colon[0,1]\to\mathbb{R}$ be almost everywhere differentiable with irrational derivative. Then, the sequence $nf(j/n)$, $1\le j\le n$ is equidistributed mod 1 in the limit $n\to\infty$.
The condition in the theorem means that there is a measure zero subset of $[0,1]$, outside of which $f$ is differentiable with $f^\prime$ irrational. Note that setting $f(x)=ax$ for a fixed irrational $a$ satisfies the requirements of the theorem. This gives Weyl's equidistribution theorem as a special case of the above result. However, this does not give an alternative approach to Weyl's theorem, as I will apply that result in the proof of the theorem above.
To answer the question asked here, we can take $f(x)=\sqrt{1-x^2}$ to show that $nf(j/n)=\sqrt{n^2-j^2}$ is equidistributed mod 1. This is differentiable on $[0,1)$ with $f^\prime$ strictly decreasing. Therefore, $f^\prime$ can only hit each rational number at most once, so can only be rational at countably many places. So, $f^\prime$ is irrational everywhere outside of a zero measure set, and the theorem above applies.
I will now prove the theorem above, the main ingredients being Van der Corput's inequality and the Weyl equidistribution theorem.
The idea of using Van der Corput's inequality comes from Terry Tao's blog post on equidistribution of polynomial sequences. However, the form given there is not quite strong enough for us. If $1\le H\le N$ are postive numbers and $\xi_n$ is a sequence of complex numbers equal to zero everywhere except $1\le n\le N$, then the following inequality holds $$ \left\lvert N^{-1}\sum_n\xi_n\right\rvert^2\le 2H^{-2}\sum_{\lvert h\rvert < H}(H-\lvert h\rvert)N^{-1}\sum_n\xi_{n+h}\bar\xi_n. $$ This is just a special case of Lemma 2.5 of Van Der Corput's Method of Exponential Sums by S.W. Graham and G. Kolesnik (link to preview containing the statement and proof of this). Given $f\colon[0,1]\to\mathbb{R}$ then, for any positive $N$, let us take $\xi_{n,N}=\exp(2\pi iNf(n/N))$ for $1\le n\le N$ and $0$ elsewhere. For any $h\ge1$, if we define $g_N\colon\mathbb{R}\to\mathbb{C}$ by $g_N(x)=\xi_{n+h}\bar\xi_n$ where $n$ is the integer satisfying $n/N\le x < (n+1)/N$ then, $$ N^{-1}\sum_{n=1}^N\xi_{n+h}\bar\xi_n=\int_{1/N}^{1+1/N}g_N(x)\,dx. $$ For any $x$ in the range $(0,1)$ then, for all $N\ge1/x$, $g(x)=\exp(2\pi i h\epsilon^{-1}(f(n/N+\epsilon)-f(n/N)))$ where $n/N\le x < n/N+\epsilon$ and $\epsilon=h/N$. So, $g_N(x)\to\exp(2\pi iNf^\prime(x))$ wherever $f$ is differentiable. If $f$ is differentiable almost everywhere then, applying bounded convergence, $$ N^{-1}\sum_{n=1}^N\xi_{n+h}\bar\xi_h\to\int_0^1\exp\left(2\pi if^\prime(x)\right)\,dx $$ as $N\to\infty$. For $h=0$ this limit is trivial, and a similar argument can be applied for $h\le-1$ (or, apply the limit above with $f(x)$ replaced by $f(1-x)$, $\xi_{n,N}$ replaced by $\xi_{N+1-n,N}$, and $h$ replaced by $-h$). So, the limit above applies for all integer $h$. Now, substituting the limit above into the Van der Corput inequality, $$ \begin{align} \limsup_{N\to\infty}\left\lvert N^{-1}\sum_n\xi_{n,N}\right\rvert^2&\le2\int_0^1\sum_{\lvert h\rvert < H}H^{-2}(H-\lvert h\rvert)\exp\left(2\pi ihf^\prime(x)\right)\,dx\\ &=2\int_0^1\sum_{h=1}^HH^{-2}\sum_{\lvert k\rvert < h}\exp\left(2\pi i kf^\prime(x)\right)\,dx. \end{align} $$ Now, for any fixed $x$ such that $f^\prime(x)$ is irrational, Weyl's equidistribution theorem says that $\sum_{\lvert k\rvert < h}\exp(2\pi ikf^\prime(x))$ is of size $o(h)$. Summing over $1\le h\le H$ then gives a term of size $o(H^2)$. So, the integrand on the right hand side of the above inequality tends to zero as $H\to\infty$ whenever $f^\prime(x)$ is irrational. If $f^\prime(x)$ is irrational almost everywhere then bounded convergence implies that the integral tends to zero. So, we have shown that $$ N^{-1}\sum_{n=1}^N\exp\left(2\pi iNf(n/N)\right)\to 0 $$ as $N\to\infty$. Applying the same result to $af(x)$ for nonzero integer $a$ shows that $$ N^{-1}\sum_{n=1}^NF(Nf(n/N))\to0=\int_0^1F(x)\,dx $$ for $F(x)=\exp(2\pi iax)$. Uniform approximation by Fourier series extends this to all continuous $F\colon\mathbb{R}\to\mathbb{C}$, showing equidistribution mod 1 of $Nf(n/N)$ over $1\le n\le N$ as $N\to\infty$.
I would state that the limit exists and it is $$\frac{1}{4}$$ I used GMP tool to see if certain hypothesis not way off the conjectured results. Numerical verification in this particular case is not a "proof" but just a navigation to the certain behavior.
I verify conjecture with GMP with 2048 bit precision floating point arithmetic. For 10k iterations it gives ~0.2479; For 20k - ~0.2487; For 40k - 0.2489; For 100k - ~0.2494
Here is a proof sketch:
Main idea here is: every continuous (and "almost everywhere" continuously differentiable - i.e. it may not be continuously differentiable in bounded number of points) function in sufficiently small area can be approximated by linear function with pre-defined precision, i.e. in $\delta, \epsilon $ - notation. This is only I have to say to state equivalence of equidistributions between current function and linear approximation. (Sorry, I do not have reference to theorem here). I.e, continuity and "almost everywhere" continuous differentiability makes equivalence.
Therefore we have to calculate limit for linear functions only.
Lets work in the first row of $\frac{1}{n}$ - type grid. Consider line going from the beginning of coordinates: $y=\frac{x}{Q}$
If Q is an even number, sum of minimum distances (between $y = 0$ and $y=\frac{1}{n}$-axes) at grid's crossing vertical lines ($x = \frac{1}{n}$) is $\frac{1}{4}$ - and it does not depends on Q.
If Q is odd number, sum is $\frac{Q^2-1}{4Q^2}$ (expression $ \to \frac{1}{4}$ when $\\Q \to +\infty$).
Let $y = \frac{P}{Q}x$, where P and Q are relatively prime. Without loss of generality let $\frac{P}{Q} < 1$. Now we can see that given set of minimum distances is just a permutation of the set of distances when P = 1 and sum depends only on Q and equal to $\frac{1}{4}$ when $\\Q \to +\infty$. To see why, lets i < Q would be the index of vertical grid line. Then we can always find $m: iP - mQ < Q$, i.e. multiplier $P$ plays role of rotational coefficient in $mod(Q)$ shift.
Now, when $y = Rx$ and R is irrational < 1, R can be approximated with rational number with arbitrary large denominator so that sum of distances $\to \frac{1}{4}$ when grid gets smaller.
That is it.
I agree that level of rigor may not be very deep, but I do not see principal problem with it in the proof.