Prove lower bound $\sum\limits_{k=1}^{n}\frac{1}{\sqrt{n^2+k^2}}\ge\left(1-\frac{1}{n}\right)\ln{(1+\sqrt{2})}+\frac{\sqrt{2}}{2n}$ 
Consider $$S_n=\sum_{k=1}^{n}\frac{1}{\sqrt{n^2+k^2}}$$ and show that, for every positive integer $n$,
  $$S_n\ge\left(1-\frac{1}{n}\right)\ln{(1+\sqrt{2})}+\dfrac{\sqrt{2}}{2n}$$

I can prove a related upper bound:
$$S_n=\sum_{k=1}^{n}\dfrac{1}{n\sqrt{1+\left(\dfrac{k}{n}\right)^2}}\le\int_{0}^{1}\dfrac{1}{\sqrt{1+x^2}}dx=\ln{(1+\sqrt{2})}$$
but I can't prove the lower bound.
 A: The inequality to be proved is 
$$\sum\limits_{k=1}^{n}\frac{1}{\sqrt{n^2+k^2}}>\left(1-\frac{1}{n}\right)\ln{(1+\sqrt{2})}+\frac{\sqrt{2}}{2n}$$
equivalent to
$$\frac{1}{n-1}\sum\limits_{k=1}^{n-1}\frac{1}{\sqrt{1+(\frac{k}{n})^2}}>\ln{(1+\sqrt{2})}=\int_0^1\frac{1}{\sqrt{1+x^2}}.$$
LHS is the Riemann Sum of integral $\displaystyle\int_0^1\frac{1}{\sqrt{1+x^2}}$. When $n\to\infty$, LHS=RHS.
So what we have to prove is that LHS is monotonically decreasing with $n$.
Consider a funtion $f(x)=\frac{1}{\sqrt{1+x^2}}+\frac{1}{\sqrt{1+(1-x)^2}}$.
Easy to know its 2nd derivative in $[0,1]$ is negative. 
So we have
$$f\Big(\frac{k}{n}\Big) > \Big(1-{k\over n}\Big)f\Big({k\over n+1}\Big)+{k\over n}f\Big({k+1\over n+1}\Big).$$ 
Sum over $n$, we get
$$\sum_{k=1}^{n-1} f\Big(\frac{k}{n}\Big)>\sum_{k=1}^{n-1}\Big(1-\frac{k}{n}\Big)f\Big(\frac{k}{n+1}\Big)+\sum_{k=1}^{n-1}\frac{k}{n}f\Big(\frac{k+1}{n+1}\Big)=\frac{n-1}{n}\sum_{k=1}^nf\Big(\frac{k}{n+1}\Big).$$
So $$\frac{1}{n-1}\sum_{k=1}^{n-1}\frac{1}{\sqrt{1+\Big(\dfrac{k}{n}\Big)^2}}=\frac{1}{2(n-1)}\sum_{k=1}^{n-1} f\Big(\frac{k}{n}\Big)$$ is monotonically decreasing with $n$. Q.E.D.
A: The bounds look like the first terms of an asymptotic series in 1/n.  I'll develop such an asymptotic series. However, this is not a complete answer to the problem, as will be stated near the end.
First, use Gradshteyn 6.611, an integral involving a Bessel function,
$ \int_0^\infty e^{-ax} J_0(bx) dx = 1/\sqrt{a^2+b^2}$ in the sum to get
$$S:=\sum_{k=1}^{n} 1/\sqrt{k^2+n^2} = \sum_{k=1}^{n}  \int_0^\infty e^{-kx} J_0(nx) dx = \int_0^\infty J_0(nx) \dfrac{1-e^{-nx}}{(e^x-1)} dx.$$
Get a series in 1/n by the following, where $B_k$ are Bernoulli numbers,
$$S=\dfrac{1}{n}  \int_0^\infty J_0(x)\dfrac{1-e^{-x}}{(e^{x/n}-1)} dx =
\int_0^\infty dx \ J_0(x)\dfrac{1-e^{-x}}{x} \sum_{k=0}^\infty \Bigl(\dfrac{x}{n}\Bigr)^k \ \dfrac{B_k}{k!}.
$$
Now do what an asymptoticist does, and interchange the sum and the integral.  The first of the two resultant integrals are actually convergent:
$$S \sim  \sum_{k=0}^\infty \Bigl(\dfrac{1}{n}\Bigr)^k \int_0^\infty J_0(x)(1-e^{-x})x^{k-1} \ dx =  $$
$$\int_0^\infty J_0(x)\frac{1-e^{-x}}{x}dx - \frac{1}{2n} \int_0^\infty J_0(x)(1-e^{-x}) + \sum_{k=1}^\infty \dfrac{B_{2k}\ n^{-2k}}{(2k)!}
\int_0^\infty dxJ_0(x)(1-e^{-x})x^{2k-1}
 $$
The first 2 integrals are $\log(1+\sqrt{2})$ and $1-1/\sqrt{2}$, respectively.
The subsequent integrals are divergent, and a method of regularization is needed for them to make sense.  I believe Chapter 4 of R. Wong, 'Asymptotic Approximation of Integrals' can be made to apply.  In particular Lemma 4 of that chapter states that 
$$\lim_{\epsilon\to0}\int_0^\infty J_\alpha(x)e^{-\epsilon \ x}x^{\mu-1}dx =
\dfrac{\Gamma(\alpha/2 + \mu/2)2^{\mu-1}}{\Gamma(\alpha/2 - \mu/2 +1)} . $$
Thus the divergent part of the integral is set to zero by a property of the $\Gamma$ function.  The portion with the exponential decay converges and can be expressed in terms of Gauss's hypergeometric function, and on collecting, we have the asymptotic series
$$S \sim \log(1+\sqrt{2}) -\frac{1-1/\sqrt{2}}{2n} - \sum_{k=1}^\infty \dfrac{B_{2k}\ n^{-2k}}{2k} F(k,k+1/2,1,-1) .$$
The standard process of summing to least term illustrates that this is a good approximation, even for small n.  For example, including only the first asymptotic term in the infinite series gives you almost five digits precision for n as small as 2!
Now for purposes of satisfying the proposer's bound question, an explicit bound on the asymptotic series must be obtained.  I have not done this, but the aforementioned chapter of Wong provides a way to do it.  The first term past those given explicitly is $-1/(24\sqrt{2} \ n^2).$  The next term of $O(n^{-4})$ is also negative.  The pattern appears to be $--++--++-- ...$ With a sufficiently precise bound on either of those first 2 negative terms, then one has proved that, for all sufficiently large n, 
$$S >   \log(1+\sqrt{2}) -\frac{1-1/\sqrt{2}}{2n} \approx 0.8814 - \dfrac{0.2929}{2n}.$$
The proposer's inequality can be written as 
$$S_p >   \log(1+\sqrt{2}) -\frac{2\log(1+\sqrt{2})-\sqrt{2}}{2n} \approx 0.8814 - \dfrac{0.3485}{2n}.$$
$S_p$ asserts a more negative $O(1/n)$ term than in $S$, so it looks like some of the negativity of the $O(1/n^2)$ asymptotic term is being pushed into the  $O(1/n)$ term.  One can probably get to $S_p$ numerically by this approach, but the closed form for the coefficient suggests a different approach. (Unless one is being sneaky, found a numerical coefficient, and then found a closed form expression that approximates the coefficient.)
A: In a previous version of this answer, we first noted that the argument to get an upper bound of $$S_n=\frac1n\sum_{k=1}^nf\left(\frac{k}n\right)\qquad f(x)=\frac1{\sqrt{1+x^2}}$$ is simply that the function $f$ is decreasing on $x>0$ hence, for every $k$, $$\frac1nf\left(\frac{k}n\right)\leqslant\int_{(k-1)/n}^{k/n}f(x)dx$$ and, summing over $1\leqslant k\leqslant n$, one gets readily $$S_n\leqslant\sum_{k=1}^n\int_{(k-1)/n}^{k/n}f(x)dx=\int_0^1f(x)dx=\left.\arg\!\sinh x\,\right|_0^1=\arg\!\sinh1=\ln(1+\sqrt2)
$$
Likewise, to get a lower bound, the most natural idea (in any case this is the idea we followed in the previous version of this answer...) might be to use once again the fact that the function $f$ is decreasing but this time, to the effect that, for every $k$, $$\frac1nf\left(\frac{k}n\right)\geqslant\int^{(k+1)/n}_{k/n}f(x)dx$$ hence, summing over $1\leqslant k\leqslant n-1$, one gets $$S_n-\frac1{\sqrt{2n^2}}\geqslant\sum_{k=1}^{n-1}\int^{(k+1)/n}_{k/n}f(x)dx=\int_{1/n}^1f(x)dx=\left.\arg\!\sinh x\,\right|_{1/n}^1=\arg\!\sinh1-\arg\!\sinh\left(\frac1n\right)
$$
Unfortunately, $$\arg\!\sinh\left(\frac1n\right)>\frac1n\arg\!\sinh1$$ hence, to get the desired lower bound, this approach is doomed and another argument is required. (The trapezoidal rule might be all that is needed here, but I did not check.)
Nota: Thanks to user @Hans for having spotted this mistake (nearly 4 years later...).
