Find triple summation rel in a closed form $S=\sum_{n=1}^{\infty}\sum_{m=1}^{n}\sum_{k=1}^{m}\frac{1}{(n+1)(k+1)(m+1)nmk}$ 
Evaluate $\displaystyle S=\sum_{n=1}^{\infty}\sum_{m=1}^{n}\sum_{k=1}^{m}\frac{1}{(n+1)(k+1)(m+1)nmk}$ 

My attempt : 
Let  $$A=\sum_{k=1}^{m}\frac{1}{k(k+1)} 
=\sum_{k=1}^{m}\left( \frac1{k}-\frac1{k+1} \right) = \frac{m}{m+1}$$
and a second sum : 
$$B=\sum_{m=1}^{n}\frac{1}{(m+1)^{2}}$$
from here how I can complete ??
 A: Continuing from where you left off, as suggested by Zarrax,
$$
\begin{align}
S&=\sum_{n=1}^\infty{1\over n(n+1)}\sum_{m=1}^n{1\over(m+1)^2}\\
&=\sum_{m=1}^\infty{1\over(m+1)^2}\sum_{n=m}^\infty{1\over n(n+1)}\\
&=\sum_{m=1}^\infty{1\over m(m+1)^2}\\
&=\sum_{m=1}^\infty\left({1\over m}-{1\over m+1}-{1\over(m+1)^2}\right)\\
&=\sum_{m=1}^\infty\left({1\over m}-{1\over m+1}\right)-\sum_{m=1}^\infty{1\over(m+1)^2}\\
&=1-\left({\pi^2\over6} -1\right)\\
&=2-{\pi^2\over6}
\end{align}
$$
A: A completely alternative approach can be *summed* up by writing the entire thing as a "single" sum:
$$S=\sum_{1\le k\le m\le n}f(k)f(m)f(n);\quad f(j)=\frac1{j(j+1)}$$
Notably, the summand is symmetric in $k,m,n$ i.e. swapping any two variables does not change the result. Using this, we can rewrite the summand with the bounds $(1\le m\le k\le n)$ or $(1\le n\le k\le m)$ etc. There are $3!=6$ ways to rearrange the bounds, but they all result in $S$. Adding every possible combination together thus gives us
$$6S=\sum_{k\ge1}\sum_{m\ge1}\sum_{n\ge1}f(k)f(m)f(n)+3\sum_{a\ge1}\sum_{b\ge1}f(a)f(b)^2+2\sum_{j\ge1}f(j)^3$$
The first triple series results from covering every possible $(k,m,n)$, since each $(k,m,n)$ must satisfy $(k\le m\le n)$ or some rearrangement of this. Then we have an extra double series for when two of the variables are equal, which has 6 possibilities, 3 of which are covered in the triple series. The last series comes about when all three variables are equal, which has 6 possibilities, 4 of which are given by the triple and double series.
Since all the indices are independent now, this is just a matter of doing some partial fractions and working each out. We have, from your work:
$$\sum_{n\ge1}f(n)=1$$
For the squared term, we have:
\begin{align}\sum_{n\ge1}f(n)^2&=\sum_{n\ge1}-\frac2n+\frac2{n+1}+\frac1{n^2}+\frac1{(n+1)^2}\\&=-\frac21-\frac1{1^2}+\sum_{n\ge1}\frac2{n^2}\\&=-3+\frac{\pi^2}3\end{align}
For the cubed term, we have:
\begin{align}\sum_{n\ge1}f(n)^3&=\sum_{n\ge1}\frac6n-\frac6{n+1}+\frac1{n^3}-\frac1{(n+1)^3}-\frac3{n^2}-\frac3{(n+1)^2}\\&=\frac61+\frac1{1^3}+\frac3{1^2}-\sum_{n\ge1}\frac6{n^2}\\&=10-\pi^2\end{align}
which finally results in
$$6S=1^3+3(1)\left(-3+\frac{\pi^2}3\right)+2(10-\pi^2)=12-\pi^2$$

$$S=2-\frac{\pi^2}6$$


More explicit explanation
Let us define
$$S_1=\sum_{k,m,n\ge1}f(k)f(m)f(n)$$
$$S_2=\sum_{k,m\ge1}f(k)f(m)^2$$
$$S_3=\sum_{k\ge1}f(k)^3$$
It should be pretty clear that $f(k)f(m)f(n)$ should occur 6 times in our $6S$. We will tackle this case-by-case.
Consider the case of $k<m<n$. As an example, consider $f(1)f(2)f(3)$. How many times does this appear in $S_1$? Notice that it occurs for
$$(k,m,n)\in\{(1,2,3),(1,3,2),(2,1,3),(2,3,1),(3,1,2),(3,2,1)\}$$
and every combination occurs in $S_1$.
Consider the case of $k\ne m=n$. As an example, consider $f(1)f(2)f(2)=f(1)f(2)^2$. How many times does this appear in $S_1$? Notice that it occurs for
$$(k,m,n)\in\{(1,2,2),(2,1,2),(2,2,1)\}$$
which is missing $3$ cases. The other $3$ case comes from $S_2$.
Consider the case of $k=m=n$. As an example, consider $f(1)f(1)f(1)=f(1)^3$. How many times does this appear in $S_1$? Notice it occurs only once, for $(k,m,n)=(1,1,1)$. How many times does this appear in $S_2$? It also occurs only once, for $(k,m)=(1,1)$. We multiply the $S_2$ by 3, which gives us 4 cases covered in total. We fill this in with 2 copies of $S_3$. This finally leaves us with
$$S=S_1+3S_2+2S_3$$
A: Yet another way:
$$\sum_{1\leq k<m<n}\frac{1}{kmn(k+1)(m+1)(n+1)} = [x^3]\prod_{h\geq 1}\left(1+\frac{x}{h(h+1)}\right)=\frac{1}{\pi}[x^4]\cos\left(\frac{\pi}{2}\sqrt{1-4x}\right)$$
where the first equality follows from Viète's theorem and the second one from the Weierstrass product for the cosine function. The Maclaurin series of the RHS immediately leads to
$$\sum_{1\leq k<m<n}\frac{1}{kmn(k+1)(m+1)(n+1)} = 5-\frac{\pi^2}{2} $$
and we may deal with the cases $k=m$ or $m=n$ in a similar fashion.
A: Continuing from the calculation above we have to sum:
$$S=\sum_{n=1}^\infty\frac{1}{n(n+1)}\sum_{m=1}^n\frac{1}{(m+1)^2}=\sum_{n=1}^{\infty}\frac{H_2(n+1)-1}{n(n+1)}\\=\sum_{n=1}^{\infty}[\frac{H_2(n+1)}{n}-\frac{H_2(n+1)}{n+1}]-\sum_{n=1}^{\infty}\frac{1}{n(n+1)}$$
Using the identity $H_2(n+1)=H_2(n)+\frac{1}{(n+1)^2}$ 
$$S=\sum_{n=1}^{\infty}[\frac{H_2(n)}{n}-\frac{H_2(n+1)}{n+1}]-1+\sum_{n=1}^{\infty}\frac{1}{n(n+1)^2}$$
And since the first sum telescopes to $1$ and the identity $\frac{1}{n(n+1)^2}=\frac{1}{n(n+1)}-\frac{1}{(n+1)^2}$ holds we obtain that:
$$S=\sum_{n=1}^{\infty}\frac{1}{n(n+1)^2}=\sum_{n=1}^{\infty}\Big[\frac{1}{n(n+1)}-\frac{1}{(n+1)^2}\Big]=1-(\zeta(2)-1)=2-\frac{\pi^2}{6}$$
A: Just for your curiosity.
We could also have a cloased form for the partial sum
$$S_p=\sum_{n=1}^{p}\sum_{m=1}^{n}\sum_{k=1}^{m}\frac{1}{(n+1)(k+1)(m+1)nmk}$$since
$$\sum_{k=1}^{m}\frac{1}{(n+1)(k+1)(m+1)nmk}=\frac{1}{(m+1)^2 n (n+1)}$$
$$\sum_{m=1}^{n}\frac{1}{(m+1)^2 n (n+1)}=\frac{\pi ^2-6-6 \psi ^{(1)}(n+2)}{6 n (n+1)}$$
$$S_p=\sum_{n=1}^{p}\frac{\pi ^2-6-6 \psi ^{(1)}(n+2)}{6 n (n+1)}=\frac{12 (p+1)-\pi ^2 (p+2)+6 (p+2) \psi ^{(1)}(p+2)}{6 (p+1)}$$
Now, using, for large values of $q$, the expansion
$$\psi ^{(1)}(q)=\frac{1}{q}+\frac{1}{2 q^2}+\frac{1}{6 q^3}+O\left(\frac{1}{q^5}\right)$$ and continuing with Taylor series for large values of $p$
$$S_p=\left(2-\frac{\pi ^2}{6}\right)+\frac{1-\frac{\pi ^2}{6}}{p}+\frac{\pi ^2-3}{6
   p^2}-\frac{2+\pi ^2}{6 p^3}+\frac{10+\pi ^2}{6 p^4}+O\left(\frac{1}{p^5}\right)$$
For illustration purposes, using $p=10$, the exact value is
$$S_{10}=\frac{42308191}{140873040}\approx 0.300328$$ while the above truncated expansion gives
$$\frac{125690-10909 \pi ^2}{60000}\approx    0.300375$$
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}$
$\ds{\mbox{Evaluate}\
\bbox[15px,#ffd]{S \equiv \sum_{n = 1}^{\infty}\sum_{m = 1}^{n}
\sum_{k = 1}^{m}{1 \over \pars{n + 1}\pars{k + 1}\pars{m + 1}nkm}}}$

\begin{align}
S & \equiv \sum_{n = 1}^{\infty}\sum_{m = 1}^{n}
\sum_{k = 1}^{m}{1 \over \pars{n + 1}\pars{k + 1}\pars{m + 1}nkm}
\\ & =
\sum_{n = 1}^{\infty}{1 \over n\pars{n + 1}}\sum_{m = 1}^{n}
{1 \over m\pars{m +1}}\
\overbrace{\sum_{k = 1}^{m}{1 \over k\pars{k + 1}}}^{\ds{m \over m + 1}}
\\ &=
\sum_{n = 1}^{\infty}{1 \over n\pars{n + 1}}\sum_{m = 1}^{n}
{1 \over \pars{m + 1}^{2}} =
\sum_{m = 1}^{\infty}{1 \over \pars{m + 1}^{2}}\
\overbrace{\sum_{n = m}^{\infty}{1 \over n\pars{n + 1}}}^{\ds{1 \over m}}
\\[5mm] & =
\sum_{m = 0}^{\infty}{1 \over \pars{m + 1}\pars{m + 2}^{2}} =
\left.-\,\totald{}{\alpha}
\sum_{m = 0}^{\infty}{1 \over \pars{m + 1}\pars{m + \alpha}}
\right\vert_{\ \alpha\ =\ 2}
\\[5mm] & =
-\,\totald{}{\alpha}\bracks{\Psi\pars{\alpha} - \Psi\pars{1} \over
\alpha - 1}
_{\ \alpha\ =\ 2}
\end{align}
$\ds{\Psi}$ is the Digamma Function and
$\ds{\Psi\pars{1} = -\gamma}$ where $\ds{\gamma}$ is the
Euler-Mascheroni Constant
Then,
\begin{align}
S & \equiv \sum_{n = 1}^{\infty}\sum_{m = 1}^{n}
\sum_{k = 1}^{m}{1 \over \pars{n + 1}\pars{k + 1}\pars{m + 1}nkm}
\\[5mm] & =
\underbrace{\Psi\pars{2}}_{\ds{1 - \gamma}}\ +\ \gamma -
\underbrace{\Psi\, '\pars{2}}_{\ds{{\pi^{2} \over 6} - 1}} =
\bbox[15px,#ffd,border:1px solid navy]{2 - {\pi^{2} \over 6}}\
\approx 0.3551
\end{align}

Note that $\ds{\Psi\pars{2}}$ is evaluated with its
recursive property $\ds{\bf\color{black}{6.3.5}}$. Namely, $\ds{\Psi\pars{2} =
\Psi\pars{\color{red}{1} + 1} = \Psi\pars{\color{red}{1}} +
{1 \over \color{red}{1}} = -\gamma + 1}$.
In addition ( see $\ds{\ \bf\color{black}{6.3.16}}$ ):
$\ds{\Psi\, '\pars{2} =
\sum_{n = 1}^{\infty}{1 \over \pars{n + 1}^{2}} =
\sum_{n = 2}^{\infty}{1 \over n^{2}} =
\sum_{n = 1}^{\infty}{1 \over n^{2}} - {1 \over 1^{2}} =
{\pi^{2} \over 6} - 1}$.
