Here's another approach.
It depends primarily on the properties of telescoping series, partial fraction expansion, and the following identity for the $m$th harmonic number
$$\begin{eqnarray*}
\sum_{k=1}^\infty \frac{1}{k(k+m)}
&=& \frac{1}{m}\sum_{k=1}^\infty \left(\frac{1}{k} - \frac{1}{k+m}\right) \\
&=& \frac{1}{m}\sum_{k=1}^m \frac{1}{k} \\
&=& \frac{H_m}{m},
\end{eqnarray*}$$
where $m=1,2,\ldots$.
Then,
$$\begin{eqnarray*}
\sum_{k=1}^{\infty}\sum_{n=1}^{\infty} \frac{1}{k^2n+2nk+n^2k}
&=& \sum_{k=1}^{\infty} \frac{1}{k}
\sum_{n=1}^{\infty} \frac{1}{n(n+k+2)} \\
&=& \sum_{k=1}^{\infty} \frac{1}{k} \frac{H_{k+2}}{k+2} \\
&=& \frac{1}{2} \sum_{k=1}^{\infty}
\left( \frac{H_{k+2}}{k} - \frac{H_{k+2}}{k+2} \right) \\
&=& \frac{1}{2} \sum_{k=1}^{\infty}
\left( \frac{H_k +\frac{1}{k+1}+\frac{1}{k+2}}{k} - \frac{H_{k+2}}{k+2} \right) \\
&=& \frac{1}{2} \sum_{k=1}^{\infty}
\left( \frac{H_k}{k} - \frac{H_{k+2}}{k+2} \right)
+ \frac{1}{2} \sum_{k=1}^{\infty}
\left(\frac{1}{k(k+1)} + \frac{1}{k(k+2)}\right) \\
&=& \frac{1}{2}\left(H_1 + \frac{H_2}{2}\right)
+ \frac{1}{2}\left(H_1 + \frac{H_2}{2}\right) \\
&=& \frac{7}{4}.
\end{eqnarray*}$$