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}$
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.
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.)