Double series convergent to $2\zeta(4)$? Using a computer I found the double sum
$$S(n)= \sum_{j=1}^n\sum_{k=1}^n \frac{j^2 + jk + k^2}{j^2(j+k)^2k^2}$$
has values
$$S(10) \quad\quad= 1.881427206538142 \\ S(1000) \quad= 2.161366028875634 \\S(100000) =  2.164613524212465\\$$
As a guess I compared with fractions $\pi^p/q$ where $p,q$ are positive integers and it appears 
$$\lim_{n \to \infty} S(n) = \frac{\pi^4}{45} =  2\zeta(4) \approx 2.164646467422276 $$
I'd be interested in seeing a proof if true.  
 A: So before I start, I've never even attempted to evaluate a double sum before, so there could very well have been an easier way.
$$\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{j^2+jk+k^2}{j^2k^2(j+k)^2} =  \sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{k^2(j+k)^2} +\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{jk(j+k)^2} + \sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2(j+k)^2}= $$
$$ 2\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{k^2(j+k)^2} +\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{jk(j+k)^2}$$
Through partial fraction decomposition the above equals:
$$\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^3k} -\frac{1}{j^3(j+k)}-\frac{1}{j^2(j+k)}+2\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{2}{j^3(j+k)}-\frac{2}{j^3k}+\frac{1}{j^2k^2}+\frac{1}{j^2(j+k)^2} $$
Collecting like-terms:
$$3\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^3(j+k)}-3\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^3k} +\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2(j+k)^2} + 2\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2k^2} =   $$
The final sum clearly equals $2\zeta(2)^2$ or $\frac{\pi^4}{18}$. I then evaluate the first two sums by combining them to get:
$$3\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^3}(\frac{1}{j+k}-\frac{1}{k}) $$
Interchanging j and k to and averaging the two sums to get:
$$\frac{3}{2}\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^3}(\frac{1}{j+k}-\frac{1}{k})+\frac{1}{k^3}(\frac{1}{j+k}-\frac{1}{j}) $$
This can be rewritten as:
$$-\frac{3}{2}\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2k(j+k)}+\frac{1}{k^2j(j+k)}= -\frac{3}{2}\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2k^2} = -\frac{3}{2}\zeta{(2)}^2 = -\frac{\pi^4}{24}$$
So putting that back into the original problem:
$$\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2(j+k)^2} + \frac{\pi^4}{18} -\frac{\pi^4}{48} = \sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2(j+k)^2} + \frac{\pi^4}{72}  $$
This is all I got to. I couldn't evaluate the other sum the way I did before. Using a calculator there is a very good chance it equals $\frac{\pi^4}{120}$.
Just for fun I was able to write the remaining sum as:
$$\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2(j+k)^2} = \sum_{n=1}^\infty \frac{\zeta(2,n+1)}{n^2}$$
Where $\zeta(x,y)$ is the Hurwitz Zeta Function. Wolfram Alpha was able to calculate the sum as $\frac{\pi^4}{120}$ as desired.
A: $$S(\infty)=\sum_{j=1}^\infty\,\sum_{k=1}^\infty \frac{(j+k)^2 - jk}{j^2(j+k)^2k^2} = \underbrace{\Big(\sum_{k=1}^\infty \frac{1}{k^2}\Big)^2}_{=\zeta(2)^2} - 
\underbrace{\sum_{j=1}^\infty\,\sum_{k=1}^\infty \frac{1}{j\,k}\int_0^\infty dt \,t \,e^{-t(j+k)}}_{:=U},$$
where the first step is algebra and the second is use of the Euler representation of the $\Gamma$ function. Interchange sums and integral and sum in terms of $\log$ to find
$$U=\int_0^\infty dt \,t \, \log^2(1-e^{-t}) = -\int_0^1 \frac{du}{u} \log\,u \log^2{(1-u)} =$$
$$-\frac{\partial}{\partial s} \frac{\partial^2}{\partial v^2} \int_0^1 u^{s-1} (1-u)^{v-1} \, du \Big\vert_{s=0,v=1}= -\frac{\partial}{\partial s} \frac{\partial^2}{\partial v^2}\frac{\Gamma(s) \Gamma(v)}{\Gamma(s+v)}\Big\vert_{s=0,v=1}$$
where the first step follows from a simple substitution $u=e^{-t}$ and the second is writing the integral in terms of something that is known, the beta integral.  Use your favorite CAS to do the partial derivatives to get $U=\pi^4/180.$ Combine with $\zeta(2)^2 = \pi^4/36$ to finish the proof of the OPs hypothesis.
A: $$ 
\begin{align} 
\frac{m^2+m\,n+n^2}{m^2\,(m+n)^2\,n^2}\, &=\frac{m^2+m\,n+n^2\,\color{red}{+m\,n-m\,n}}{m^2\,(m+n)^2\,n^2} \\[2mm] 
&=\,\frac{(m+n)^2-m\,n}{(m+n)^2\,m^2\,n^2} \\[2mm] 
&=\,\frac{1}{m^2\,n^2}-\frac{1}{m\,n\,(m+n)^2} \\[2mm] 
&=\,\frac{1}{m^2\,n^2}-\frac{1}{m^3}\left(\frac{1}{n}-\frac{1}{m+n}-\frac{m}{(m+n)^2}\right) \\[2mm] 
&=\,\color{brown}{\frac{1}{m^2\,n^2}}\color{green}{-\frac{1}{m^3}\left(\frac{1}{n}-\frac{1}{m+n}\right)}\color{blue}{+\frac{1}{m^2}\frac{1}{(m+n)^2}} 
\end{align} 
$$ 

$$ 
\begin{align} 
\color{brown}{\large S_{\small 1}\,} &=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\frac{1}{m^2\,n^2}=\sum_{m=1}^{\infty}\frac{1}{m^2}\sum_{n=1}^{\infty}\frac{1}{n^2}=\left(\zeta(2)\right)^2=\color{brown}{\frac{\pi^4}{36}} \\[4mm] 
\color{green}{\large S_{\small 2}\,} &=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\frac{1}{m^3}\left(\frac{1}{n}-\frac{1}{m+n}\right)=\sum_{m=1}^{\infty}\frac{1}{m^3}\sum_{n=1}^{\infty}\left(\frac{1}{n}-\frac{1}{m+n}\right) \\[2mm] 
&=\sum_{m=1}^{\infty}\frac{1}{m^3}\sum_{n=1}^{\color{red}{m}}\frac{1}{n}=\sum_{m=1}^{\infty}\frac{H_{m}}{m^3}=\frac{5}{4}\zeta(4)=\color{green}{\frac{\pi^4}{72}}\tag{1} \\[4mm] 
\color{blue}{\large S_{\small 3}\,} &=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\frac{1}{m^2}\frac{1}{(m+n)^2}=\sum_{m=1}^{\infty}\frac{1}{m^2}\sum_{n=1}^{\infty}\frac{1}{(m+n)^2} \\[2mm] 
&=\sum_{m=1}^{\infty}\frac{1}{m^2}\,\sum_{\color{red}{n=m+1}}^{\infty}\,\frac{1}{n^2}=\sum_{m=1}^{\infty}\frac{\psi^{\small (1)}(m+1)}{m^2} \\[2mm] 
&=\sum_{m=1}^{\infty}\frac{1}{m^2}\left[\zeta(2)-\sum_{k=1}^{m}\frac{1}{k^2}\right]=\left(\zeta(2)\right)^2-\sum_{m=1}^{\infty}\frac{H_{m,2}}{m^2} \\[2mm] 
&=\left(\zeta(2)\right)^2-\frac{1}{2}\left[\left(\zeta(2)\right)^2+\zeta(4)\right]=\frac{1}{2}\left[\left(\zeta(2)\right)^2-\zeta(4)\right]=\color{blue}{\frac{\pi^4}{120}}\tag{2} 
\end{align} 
$$ 
$$ \color{red}{\Longrightarrow\quad S}\,=S_1-S_2+S_3=\,\color{red}{\frac{\pi^4}{45}} $$ 

$\,H_m\,\,\,$ : Harmonic Number , $\,\{1\}\,$ : Equation (19) ,$\,\{2\}\,$ : Equation (43) 
$\,{\large\psi}^{\small (1)}\,\,$ : Polygamma Function 
A: An alternative approach:
$$ S = \lim_{n\to +\infty}S(n) = \sum_{j,k\geq 1}\frac{1}{j^2 k^2}-\sum_{k,j\geq 1}\frac{1}{jk(j+k)^2}=\zeta(2)^2-\sum_{k,j\geq 1}\int_{0}^{+\infty}\frac{e^{-(j+k)x}}{jk}\,x\,dx $$
leads to
$$S = \zeta(2)^2-\int_{0}^{+\infty}x\log^2(1-e^{-x})\,dx=\frac{\pi^4}{36}+\int_{0}^{1}\frac{\log^2(1-x)\log(x)}{x}\,dx$$
or to
$$ S = \frac{\pi^4}{36}+\int_{0}^{1}\frac{\log(1-x)}{1-x}\log^2(x)\,dx = \frac{\pi^4}{36}-\sum_{n\geq 1}\frac{2H_n}{(n+1)^3}$$
since $\frac{-\log(1-x)}{1-x}=\sum_{n\geq 1}H_n x^n$ and $\int_{0}^{1}x^n\log^2(x)\,dx = \frac{2}{(n+1)^3}$. Rearranging
$$ S = \frac{\pi^4}{36}-2\sum_{n\geq 1}\frac{H_{n}}{n^3}+2\,\zeta(4) = 2\,\zeta(4) = \frac{\pi^4}{45}$$
since the middle term is a linear Euler sum, which can be computed from the Theorem 2.2 here (Flajolet and Salvy, a masterpiece).
A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align}
&\bbox[#ffd,10px]{\ds{\sum_{j = 1}^{\infty}\sum_{k = 1}^{\infty}
{j^{2} + jk + k^{2} \over j^{2}\pars{j + k}^{2}k^{2}}}} =
\sum_{j = 1}^{\infty}\sum_{k = 1}^{\infty}
{j^{2} + jk + k^{2} \over j^{2}k^{2}}\
\overbrace{\bracks{-\int_{0}^{1}\ln\pars{x}x^{j + k - 1}\,\dd x}}
^{\ds{1 \over \pars{j + k}^{2}}}
\\[5mm] = &\
-\int_{0}^{1}\ln\pars{x}
\pars{\sum_{j = 1}^{\infty}x^{j}
\sum_{k = 1}^{\infty}{x^{k} \over k^{2}} +
\sum_{j = 1}^{\infty}{x^{j} \over j}
\sum_{k = 1}^{\infty}{x^{k} \over k} +
\sum_{j = 1}^{\infty}{x^{j} \over j^{2}}
\sum_{k = 1}^{\infty}x^{k}}\,{\dd x \over x}
\\[5mm] = &\
-\int_{0}^{1}\ln\pars{x}
\bracks{2\sum_{j = 1}^{\infty}x^{j}
\sum_{k = 1}^{\infty}{x^{k} \over k^{2}} +
\pars{\sum_{j = 1}^{\infty}{x^{j} \over j}}^{2}}\,{\dd x \over x}
\label{1}\tag{1}
\end{align}

Note that $\ds{\sum_{\ell = 1}^{\infty}{x^{\ell} \over \ell^{s}} =
\,\mrm{Li}_{s}\pars{x}}$ where $\ds{\mrm{Li}_{s}}$ is the Polylogarithm Function. Moreover,
  $\ds{\mrm{Li}_{1}\pars{x} = -\ln\pars{1- x}}$,
  $\ds{\mrm{Li}_{s + 1}\pars{z} = \int_{0}^{z}{\mrm{Li}_{s}\pars{t} \over t}\,\dd t}$ and
  $\ds{\sum_{j = 1}^{\infty}x^{j} = {x \over 1 - x}}$.

\eqref{1} becomes
\begin{align}
&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\bbox[#ffd,10px]{\ds{\sum_{j = 1}^{\infty}\sum_{k = 1}^{\infty}
{j^{2} + jk + k^{2} \over j^{2}\pars{j + k}^{2}k^{2}}}} =
-2\int_{0}^{1}{\ln\pars{x}\,\mrm{Li}_{2}\pars{x} \over 1 - x}\,\dd x -
\int_{0}^{1}{\ln\pars{x}\ln^{2}\pars{1 - x} \over x}\,\dd x
\\[5mm] = &\
-2\int_{0}^{1}\overbrace{\ln\pars{1 - x} \over x}^{\ds{-\mrm{Li}_{2}'\pars{x}}}\,\
\mrm{Li}_{2}\pars{1 - x} \,\dd x -
\int_{0}^{1}{\ln\pars{x}\ln^{2}\pars{1 - x} \over x}\,\dd x
\label{2}\tag{2}
\end{align}
where we set $\ds{x \mapsto 1 - x}$ in the first integral. With Euler Reflection Formula
$\ds{\mrm{Li}_{2}\pars{1 - x} = -\mrm{Li}_{2}\pars{x} + {\pi^{2} \over 6} -\ln\pars{x}\ln\pars{1 - x}}$. \eqref{2} becomes:
\begin{align}
&\bbox[#ffd,10px]{\ds{\sum_{j = 1}^{\infty}\sum_{k = 1}^{\infty}
{j^{2} + jk + k^{2} \over j^{2}\pars{j + k}^{2}k^{2}}}}
\\ = &\ -\
\overbrace{\int_{0}^{1}\totald{\mrm{Li}_{2}^{2}\pars{x}}{x} \,\dd x}
^{\ds{\pi^{4} \over 36}}\ +\
{\pi^{2} \over 3}\
\overbrace{\int_{0}^{1}\mrm{Li}_{2}'\pars{x}\,\dd x}
^{\ds{\pi^{2} \over 6}}\ +\
\overbrace{\int_{0}^{1}{\ln\pars{x}\ln^{2}\pars{1 - x} \over x}\,\dd x}^{\ds{-\,{\pi^{4} \over 180}}}
\\[5mm] = &\
\bbx{\large{\pi^{4} \over 45}}
\end{align}
