How to get$$\mathop {\lim }\limits_{n \to \infty } n\left( {\frac{\pi }{2} + \ln 2 - \frac{1}{n}\sum\limits_{j = 1}^n {\sum\limits_{i = 1}^n {\frac{{i + j}}{{{i^2} + {j^2}}}} } } \right).$$
I think we can use Euler–Maclaurin formula$$\sum_{n=a}^b f(n) \sim \int_a^b f(x)\,\mathrm{d}x + \frac{f(b) + f(a)}{2} + \sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} \left(f^{(2k - 1)}(b) - f^{(2k - 1)}(a)\right),$$ where $a,b$ are both integers. But it seems difficult because of the double summation!
Actually, for $S_n=\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{i+j}{i^2+j^2}$, the limit $S=\displaystyle\lim_{n\to\infty}\left[\left(\frac{\pi}{2}+\ln 2\right)n\color{red}{-\ln n}-S_n\right]$ exists.
This limit is seen as $1+\displaystyle\lim_{n\to\infty}(I_n-S_n)$, where $$\begin{align*}I_n&=\iint_{[\frac{1}{2},n+\frac{1}{2}]^2}\frac{x+y}{x^2+y^2}\,dx\,dy\\ &=(2n+1)\ln(2n+1)\\ &-(n+1)\ln(2n^2+2n+1)\\ &+n(2\arctan(2n+1)-\pi/2).\end{align*}$$ To prove existence, note that $I_n-S_n=\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}\Delta_{i,j}$, where $$\begin{gather}\Delta_{i,j}=\iint_{[-\frac{1}{2},\frac{1}{2}]^2}\big(f(i+x,j+y)-f(i,j)\big)\,dx\,dy,\quad f(x,y)=\frac{x+y}{x^2+y^2};\\ \frac{1}{n!}\left|\frac{\partial^{n}\!f}{\partial x^k\partial y^{n-k}}\right|\leqslant\left|\frac{x+y}{x^2+y^2}\right|^{n+1},\quad \frac{\partial^2\!f}{\partial x^2}+\frac{\partial^2\!f}{\partial y^2}=0,\end{gather}$$ and Taylor's theorem produces $|\Delta_{i,j}|=\mathcal{O}\left(\big(\frac{i+j}{i^2+j^2}\big)^5\right)$ which is sufficient.
Computing $S$ using (its definition and) Lagrange-Zagier extrapolation, I get $$ \color{blue}{S=1.0042628439817233943074076864477736788445647263436\ldots} $$ I wonder whether $S$ is related to some known mathematical constants...
The higher-order asymptotics can indeed be derived from Euler-Maclaurin formula (or its two-dimensional extension, but I'm going the elementary way below). Let $$ \begin{gather}S_n\asymp Kn-\ln n-S-\sum_{k=1}^{(\infty)}\frac{a_k}{n^k},\quad K=\frac{\pi}{2}+\ln 2;\\ S_n-S_{n-1}\asymp K+\ln\Big(1-\frac{1}{n}\Big)+\sum_{k=2}^{(\infty)}n^{-k}\sum_{r=1}^{k-1}\binom{k-1}{r-1}a_r.\end{gather} $$ Now we apply E.-M. to $f(x)=2\dfrac{1+x}{1+x^2}, x\in[0,1]$. We have $\displaystyle\int_0^1 f(x)\,dx=K$, $$ \begin{gather}\frac{1}{n}\left(\frac{f(0)+f(1)}{2}+\sum_{k=1}^{n-1}f\Big(\frac{k}{n}\Big)\right)=\frac{1}{n}+S_n-S_{n-1},\\ \frac{f^{2k-1}(0)}{(2k)!}=\frac{(-1)^{k-1}}{k},\quad\quad\frac{f^{2k-1}(1)}{(2k)!}=-\frac{(-1)^{\lfloor k/2\rfloor}}{2^k\cdot k}\end{gather} $$ (the last two e.g. from power series); Euler-Maclaurin gives $$ S_n-S_{n-1}\asymp K-\frac{1}{n}-\sum_{k=1}^{(\infty)}\frac{c_k}{n^{2k}},\quad c_k=\frac{B_{2k}}{k}\Big((-1)^{k-1}+2^{-k}(-1)^{\lfloor k/2\rfloor}\Big). $$ Thus, $$ \sum_{r=1}^{k-1}\binom{k-1}{r-1}a_r=b_k=\frac{1}{k}-\begin{cases}c_{k/2}&(k\text{ even})\\ 0 &(k\text{ odd})\end{cases}\quad(k>1) $$ and, recognizing the inverse matrix for this system here, we finally get $$ a_n=\frac{1}{n}\sum_{k=1}^{n}\binom{n}{k}B_{n-k}b_{k+1}. $$ This sequence begins with $$ \frac{1}{4}, \frac{1}{24}, -\frac{7}{144}, \frac{3}{160}, 0, -\frac{1}{2016}, -\frac{19}{2688}, \frac{31}{3840}, 0, \frac{1}{4224}, -\frac{1453}{59136}, \frac{29713}{698880}, 0, \ldots $$ (this coincides with the values computed numerically in a prior edition of this answer).
The main purpose of this post is the analytical evaluation of the constant term.
Let's denote $$ S(n)=\sum_{k=1}^n\sum_{l=1}^n\frac{k+l}{k^2+l^2}=2\sum_{k=1}^n\sum_{l=1}^n\frac{l}{k^2+l^2}=2\sum_{k=1}^\infty\sum_{l=1}^n\frac{l}{k^2+l^2}-2\sum_{k=n+1}^\infty\sum_{l=1}^n\frac{l}{k^2+l^2}$$ $$=S_1+S_2$$ To evaluate $S_1$ we change the order of summation $$S_1=2\sum_{l=1}^n\sum_{k=1}^\infty\frac{l}{k^2+l^2}=\sum_{l=1}^n\sum_{k=-\infty}^\infty\frac{l}{k^2+l^2}-\sum_{l=1}^n\frac{1}{l}$$ Using $\displaystyle \sum_{k=-\infty}^\infty\frac{1}{k^2+l^2}=\frac{\pi}{l}\coth\pi l=\frac{\pi}{l}\Big(1+\frac{2}{e^{2\pi l}-1}\Big)$ $$S_1=\pi n-\sum_{l=1}^n\frac{1}{l}+\sum_{l=1}^n\frac{2\pi}{e^{2\pi l}-1}$$ With the accuracy up to exponentially small terms we can expand summation to $\infty$. Using also the asymptotics of the second term at $n\to\infty$, we can get as many terms as we want.
For our purpose, $$S_1=\pi n-\ln n-\gamma+\sum_{l=1}^n\frac{2\pi}{e^{2\pi l}-1}+O\Big(\frac{1}{n}\Big)\tag{1}$$ Now, evaluating $\,S_2$ $$S_2=-2\,\Im\sum_{k=n+1}^\infty\sum_{l=1}^n\frac{1}{k-il}=-2\,\Im S_0(n)\tag{2}$$ where $\displaystyle S_0(n)=\sum_{k=n+1}^\infty\sum_{l=1}^n\int_0^\infty e^{-v(k-il)}dv$
Changing the order of summation and integration $$S_0=\int_0^\infty \frac{e^{-v(n+1)}}{1-e^{-v}}\frac{e^{iv}-e^{iv(n+1)}}{1-e^{iv}}dv=\frac{i}{2}\int_0^\infty \frac{e^{-vn}}{e^v-1}\frac{e^{iv/2}-e^{ivn+iv/2}}{\sin\frac{v}{2}}dv$$ $$S_2=-2\,\Im S_0(n)=-\int_0^\infty\frac{e^{-vn}}{e^v-1}\Big(\cot\big(\frac{v}{2}\big)\big(1-\cos vn\big)+\sin vn\Big)dv$$ $$=-\frac{1}{n}\int_0^\infty\frac{e^{-t}}{e^{t/n}-1}\Big(\cot\big(\frac{t}{2n}\big)\big(1-\cos t\big)+\sin t\Big)dv\tag{3}$$ Due to the exponent in the numerator we can decompose the integrand into the series of powers $\frac{1}{n}$. Again, we can evaluate as many terms of the decomposition as we want.
With the desired accuracy in our case $$S_2=-2n\int_0^\infty\frac{1-\cos t}{t^2}e^{-t}dt+\int_0^\infty\frac{1-\cos t}{t}e^{-t}dt-\int_0^\infty\frac{\sin t}{t}e^{-t}dt+O\Big(\frac{1}{n}\Big)$$ Integration is straightforward: $$S_2=\Big(\ln 2-\frac{\pi}{2}\Big)n+\frac{\ln 2}{2}-\frac{\pi}{4}+O\Big(\frac{1}{n}\Big)\tag{4}$$ Taking together (1) and (4) $$\boxed{\,\,S=S_1+S_2=\Big(\frac{\pi}{2}+\ln 2\Big)n-\ln n+\bigg(\sum_{l=1}^n\frac{2\pi}{e^{2\pi l}-1}-\gamma+\frac{\ln 2}{2}-\frac{\pi}{4}\bigg)+O\Big(\frac{1}{n}\Big)\,\,}$$ and the desired constant term $\,\,\displaystyle\sum_{l=1}^\infty\frac{2\pi}{e^{2\pi l}-1}-\gamma+\frac{\ln 2}{2}-\frac{\pi}{4}=-1.00426..$
