Playing around with double sums related to harmonic numbers and having easily found a closed expression for this sum
$$s_1 = \sum_{n,m\ge 1} \frac{1}{n m (n+m)}=2 \zeta(3)\tag{1}$$
I attempted to find a closed expression for this slightly more complicated sum
$$s_2 = \sum_{n,m\ge 1} \frac{1}{n m (n^2+m^2)}\tag{2}$$
which turned out to be more than I could achieve, and I hope for your help.
Here's what I did so far.
First approach
$$\begin{align}s_2 & = \sum_{n,m\ge 1} \frac{1}{n m (n^2+m^2)} \\ & = \sum_{n,m\ge 1} \frac{1}{n m }\int_{0}^{\infty } e^{-t(n^2+m^2)}\;dt\\ & =\int_{0}^{\infty }\left(\sum_{n,m\ge 1} \frac{1}{n m } e^{-t(n^2+m^2)}\right)\;dt\\ & =\int_{0}^{\infty } f_2(t)^2\;dt \end {align}\tag{3}$$
where
$$f_2(t) = \sum_{n\ge 1}\frac{e^{-t n^2}}{n}\tag{4}$$
Now
$$e^{-n^2 t}=\frac{1}{\sqrt{\pi}}\int_{-\infty }^{\infty } e^{-x^2+i 2 n x\sqrt{t} } \, dx\tag{5}$$
This linearizes $n$ in the exponent, and, exchanging integral and sum, we can do the (conditionally convergent) $n$-sum under the $x$-integral which gives
$$\frac{1}{\sqrt{\pi}}\sum _{n=1}^{\infty } \frac{e^{-x^2+i 2 n x \sqrt{t}}}{\sqrt{\pi } n}=-\frac{e^{-x^2}}{\sqrt{\pi}} \log \left(1-e^{2 i x \sqrt{t} }\right)\tag{6}$$
and
$$f_2(t) = -\frac{1}{\sqrt{\pi}}\int_{-\infty }^{\infty } e^{-x^2} \log \left(1-e^{2 i x \sqrt{t} }\right)\, dx\tag{7}$$
Now in the $x$-integral the imaginary part cancels out due to symmetry and the real part can be simplified so that
$$f_2(t) = -\frac{1}{\sqrt{\pi}}\int_{-\infty }^{\infty } e^{-x^2} \frac{1}{2} \log \left(4 \sin ^2\left(x\sqrt{t} \right)\right)\; dx\tag{8}$$
Alas, here I am stuck, and summing up it seems that I have made a simple looking formula more complicated.
Second approach
Doing just the $n$-sum in $s_2$ gives
$$s_2 = \sum_{m\ge 1} \frac{1}{2 m^3} (H(i m) + H(-i m))=\Re\left(\sum_{m\ge 1} \frac{1}{m^3} H(i m)\right)\tag{2.1}$$
where $H(z)$ is the harmonic number of argument $z$.
Inserting the representation
$$H(z) = \int_{0}^{1} \frac{1-x^z}{1-x}\;dx$$
and performing the $m$-sum gives an interesting integral representation of our double sum
$$s_2 = \left(\int_0^1 \frac{-\zeta (3)+\Re(\text{Li}_3\left(x^i\right))}{x-1} \, dx\right)\tag{2.2}$$