Inspired by Find the sum of the double series $\sum_{k=1}^\infty \sum_{j=1}^\infty \frac{1}{(k+1)(j+1)(k+1+j)} $ here are some related problems.
- Prove that
$$w(2,1) = \sum _{i=1}^{\infty } \sum _{j=1}^{\infty } \frac{1}{i j (i+j)} = 2 \zeta (3)$$
$$w(3,1)=\sum _{i=1}^{\infty } \sum _{j=1}^{\infty } \sum _{k=1}^{\infty } \frac{1}{i j k (i+j+k)}=6 \zeta (4)$$
And can you guess (and prove) the closed form result for
$$w(n,1)=\sum _{i_{1}=1}^{\infty } \sum _{i_{2}=1}^{\infty }... \sum _{i_{n}=1}^{\infty } \frac{1}{i_{1} i_{2}... i_{n} (i_{1}+i_{2}+...+i_{n})}$$
- Calculate, if possible, closed forms for
$$w(2,2,1) = \sum _{i=1}^{\infty } \sum _{j=1}^{\infty } \frac{1}{i^2 j^2 (i+j)} $$
$$w(2,2,2) = \sum _{i=1}^{\infty } \sum _{j=1}^{\infty } \frac{1}{i^2 j^2 (i+j)^2} $$
The question of the OP was completely answered by Jack D'Aurizio. Now, as an extension, I'd like to study more general problems of a similar type.
$$w(p,q,r)=\sum _{n=1}^{\infty } \sum _{m=1}^{\infty } \frac{1}{m^p n^q (m+n)^r}\tag{1}$$
and start with the case $p = q$ general $r$ and see how far we can come.
We will consider small values in order to see a pattern.
The solution method will be direct summation, exchange of order of summation and partial fraction decomposition. I discarded use of integrals which did not lead me to results.
$$H_k=\sum _{i=1}^k \frac{1}{i}$$ $$H_k^{(m)}=\sum _{i=1}^k \frac{1}{i^m}$$ $$\zeta (s)=\sum _{n=1}^{\infty } \frac{1}{n^s}$$
$$S(1,q)/;(q\ge 2) = \sum _{k=1}^{\infty } \frac{H_k}{k^q}= \left(\frac{q}{2}+1\right) \zeta (q+1)-\frac{1}{2} \sum _{k=1}^{q-2} \zeta (k+1) \zeta (q-k)\tag{2}$$
$$S(p,q)/;(p+q\; odd)=\sum _{k=1}^{\infty } \frac{H_k^{(p)}}{k^q}=(-1)^p \sum _{k=1}^{\left\lfloor \frac{q}{2}\right\rfloor } \zeta (2 k) \binom{-2 k+p+q-1}{p-1} \zeta (-2 k+p+q)+(-1)^p \sum _{k=1}^{\left\lfloor \frac{p}{2}\right\rfloor } \zeta (2 k) \binom{-2 k+p+q-1}{q-1} \zeta (-2 k+p+q)+\left(-\frac{1}{2} (-1)^p \binom{p+q-1}{p}-\frac{1}{2} (-1)^p \binom{p+q-1}{q}+\frac{1}{2}\right) \zeta (p+q)+\frac{1}{2} \left(1-(-1)^p\right) \zeta (p) \zeta (q)\tag{3}$$
$$ \begin{eqnarray*} w(p,p,r)&=&\sum _{n=1}^{\infty } \sum _{m=1}^{\infty } \frac{1}{m^p n^p (m+n)^r}\\ &=&\sum _{n} \sum _{s}\frac{1}{(s-n)^p n^p s^r}\\ &=&\sum _{s \ge 2} \frac{1}{s^r}{\sum _{n=1}^{r-1} \frac{1}{(n (s-n))^p}}\\ &=&\sum _{n} \sum _{s}\frac{1}{(s-n)^p n^p s^r}\\ &=&\sum _{s \ge 2} \frac{1}{s^{r+p}}{\sum _{n=1}^{r-1} (\frac{1}{n}+\frac{1}{ (s-n)})^p}\tag{4a}\\ &=&\sum _{s \ge 2} \frac{1}{s^{p+r}} \sum _{n=1}^{r-1} \left(\sum _{t=0}^p \binom{p}{t} \frac{1}{n^t (s-n)^{p-t}}\right)\tag{4b}\\ \end{eqnarray*} $$
Now let us calculate the first few p-sums.
$p= 1$
$$ \begin{align} w(1,1,r) &= \sum _{s=2}^{\infty } \frac{1}{s^{r+1}}\sum _{n=1}^{r-1} \left(\frac{1}{s-n}+\frac{1}{n}\right)=2\sum _{s=2}^{\infty } \frac{ H_{s-1}}{s^{r+1}}\\ &= 2 \sum _{s=2}^{\infty } \frac{H_s-\frac{1}{s}}{s^{r+1}}=2 \sum _{s=2}^{\infty } \frac{H_s}{s^{r+1}}-2 \;\zeta(r+2)\\ &= \left(r+1\right) \zeta (r+2)- \sum _{k=1}^{r-1} \zeta (k+1) \zeta (-k+r+1)\tag{5} \end{align} $$
Hence for $p=1$ the sum reduces to zeta values for any $r\ge 2$.
For the first 5 values of r we have
$$ \begin{array}{l} w(1,1,1)=2 \zeta (3) \\ w(1,1,2)= 3 \zeta(4)-\zeta(2)^2=\frac{\pi ^4}{180} \\ w(1,1,3) = 4 \zeta(5)-2 \zeta(2) \zeta(3)=4 \zeta (5)-\frac{\pi ^2 \zeta (3)}{3} \\ w(1,1,4) = 5 \zeta(6)-\zeta(3)^2-2 \zeta(2) \zeta(4)=\frac{\pi ^6}{630}-\zeta (3)^2 \\ w(1,1,5) = 6 \zeta(7) -2 \zeta(3) \zeta(4)-2 \zeta(2) \zeta(5)=-\frac{\pi ^4 \zeta (3)}{45}-\frac{\pi ^2 \zeta (5)}{3}+6 \zeta (7) \\ \end{array} $$
$p= 2$
Since
$$\left(\frac{1}{s-n}+\frac{1}{n}\right)^2=\frac{1}{n^2}+\frac{2}{n (s-n)}+\frac{1}{(s-n)^2}=\frac{1}{n^2}+\frac{2}{s} \left(\frac{1}{s-n}+\frac{1}{n}\right)+\frac{1}{(s-n)^2}$$
the n-sum becomes
$$\sum _{n=1}^{s-1} \left(\frac{1}{n^2}+\frac{2}{s} \left(\frac{1}{s-n}+\frac{1}{n}\right)+\frac{1}{(s-n)^2}\right)=\frac{4 H_{s-1}}{s}+2 H_{s-1}^{(2)}$$
Hence
$$ \begin{align} w(2,2,r)&=4 \sum _{s=1}^{\infty } \frac{H_{s-1}}{s^{3+r}}+2 \sum _{s=1}^{\infty } \frac{H_{s-1}^{(2)}}{s^{2+r}}\\ &=-6 \zeta(4+r)+4 \sum _{s=1}^{\infty } \frac{H_{s}}{s^{3+r}}+2 \sum _{s=1}^{\infty } \frac{H_{s}^{(2)}}{s^{2+r}}\tag{6a}\\ &= -6 \zeta(4+r)+ 4 S(1,3+r)+ 2 S(2,2+r)\tag{6b} \end{align} $$
Notice that in $(6b)$ the explicit formula $(3)$ for the sums $S(p,q)$ are only valid for odd $r$.
Hence for $p=2$ and odd $r$ we have a solution in terms of zeta values.
For even $r$ there seems to be no general formula for S(2,2m) (more generally for even p+q) and we have to take $(6a)$
The "problem child" is the sum
$$\sum _{s=1}^{\infty } \frac{H_{s}^{(2)}}{s^{2+r}}/; r \; even$$
I have only seen these two cases which the even sums reduce to zeta values (formulas (40) and (42) resp. in [1])
$$S(2,4)=\zeta (3)^2-\frac{\pi ^6}{2835}$$ $$S(4,2)=\frac{37 \pi ^6}{11340}-\zeta (3)^2$$
$p=3$
to be done.
References
[1] http://mathworld.wolfram.com/HarmonicNumber.html
[2] http://algo.inria.fr/flajolet/Publications/FlSa98.pdf
[3] http://mathworld.wolfram.com/EulerSum.html