Asymptotic estimation problem about $\sum\limits_{j = 1}^n {\sum\limits_{i = 1}^n {\frac{{i + j}}{{{i^2} + {j^2}}}} } $ 
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!
 A: 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{l}{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..$
