Mysterious subleading corrections to sums with internal dependence on limit
Solution 1:
Usually one can use the Euler-Maclaurin sum formula to get information like this. Its first iteration takes the form
$$ \begin{align} &\frac{1}{n} \sum_{k=0}^{n-1} f\left(\frac{k}{n}\right) \\ &\qquad = \int_0^1 f(x)\,dx + \frac{f(0)-f(1)}{2n} - \frac{1}{n} \int_0^1 \left(-nx - \lceil -nx\rceil - \frac{1}{2}\right)f'(x)\,dx, \end{align} $$
and when applied to your sum we get
$$ \frac{S(n)}{n^2} = \frac{\pi}{4} + \frac{1}{2n} + \frac{1}{n} \int_0^1 \left(-nx - \lceil -nx\rceil - \frac{1}{2}\right)\frac{x}{\sqrt{1-x^2}}\,dx. $$
Now it just remains to estimate the integral. To do this, we first expand the periodic factor in Fourier series as
$$ -nx - \lceil -nx\rceil - \frac{1}{2} = \frac{1}{\pi} \sum_{k=1}^{\infty} \frac{\sin(2\pi k n x)}{k}, $$
so that
$$ \begin{align} &\int_0^1 \left(-nx - \lceil -nx\rceil - \frac{1}{2}\right)\frac{x}{\sqrt{1-x^2}}\,dx \\ &\qquad = \frac{1}{\pi} \sum_{k=1}^{\infty} \frac{1}{k} \int_0^1 \sin(2\pi k n x) \frac{x}{\sqrt{1-x^2}}\,dx \\ &\qquad = \frac{1}{\pi} \sum_{k=1}^{\infty} \frac{1}{k} \left[\left. -\sin(2\pi k n x) \sqrt{1-x^2} \right|_0^1 + 2\pi k n \int_0^1 \cos(2\pi k n x) \sqrt{1-x^2}\,dx\right] \\ &\qquad = \frac{1}{\pi} \sum_{k=1}^{\infty} \frac{1}{k} \cdot 2\pi k n \int_0^1 \cos(2\pi k n x) \sqrt{1-x^2}\,dx \\ &\qquad = \frac{1}{2} \sum_{k=1}^{\infty} \frac{J_1(2\pi k n)}{k}, \end{align} $$
where $J_1$ is the Bessel function of the first kind of order $1$. This Bessel function has the asymptotic form
$$ J_1(x) = - \sqrt{\frac{2}{\pi x}} \cos\left(x+\frac{\pi}{4}\right) + O\left(x^{-3/2}\right) $$
as $x \to \infty$, and since $\cos(2\pi k n + \pi/4) = 1/\sqrt{2}$ we get
$$ J_1(2\pi k n) = - \frac{1}{\pi \sqrt{2 k n}} + O\left((kn)^{-3/2}\right). $$
It follows that
$$ \begin{align} \sum_{k=1}^{\infty} \frac{J_1(2\pi k n)}{k} &= -\frac{1}{\pi \sqrt{2n}} \sum_{k=1}^{\infty} \frac{1}{k^{3/2}} + O\left(n^{-3/2}\right) \\ & = -\frac{\zeta(3/2)}{\pi\sqrt{2n}} + O\left(n^{-3/2}\right) \end{align} $$
and hence that
$$ \int_0^1 \left(\lceil nx\rceil - nx - \frac{1}{2}\right)\frac{x}{\sqrt{1-x^2}}\,dx = -\frac{\zeta(3/2)}{\pi2^{3/2}\sqrt{n}} + O\left(n^{-3/2}\right). $$
Therefore
$$ \frac{S(n)}{n^2} = \frac{\pi}{4} + \frac{1}{2n} - \frac{\zeta(3/2)}{\pi 2^{3/2} n^{3/2}} + O\left(n^{-5/2}\right). $$
or
$$ S(n) = \frac{\pi}{4}n^2 + \frac{1}{2}n - \frac{\zeta(3/2)}{\pi 2^{3/2}} n^{1/2} + O\left(n^{-1/2}\right). $$
Note that the constant you estimated on the $n^{1/2}$ term is
$$ - \frac{\zeta(3/2)}{\pi 2^{3/2}} \approx -0.293\ 995\ 518\ 793\ 519\ 291. $$