Is there a closed-form solution for $\sum_{n=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{nm(3n+m)}$? I am seeking a closed-form solution for this double sum:
\begin{eqnarray*}
\sum_{n=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{nm(\color{blue}{3}n+m)}= ?.
\end{eqnarray*}
I will turn it into $3$ tough integrals in a moment.  But first I will state some similar results:
\begin{eqnarray*}
\sum_{n=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{nm(n+m)} &=& 2 \zeta(3) \\
\sum_{n=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{nm(\color{blue}{2}n+m)} &=& \frac{11}{8} \zeta(3) \\
\sum_{n=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{nm(\color{blue}{4}n+m)} &=& \frac{67}{32} \zeta(3) -\frac{G \pi}{2}. \\
\end{eqnarray*}
where $G$ is the Catalan constant. The last result took some effort ... 
Now I know most of you folks prefer integrals to sums, so lets turn this into an integral. Using
\begin{eqnarray*}
\frac{1}{n} &=& \int_0^1 x^{n-1} dx\\
\frac{1}{m} &=& \int_0^1 y^{m-1} dy\\
\frac{1}{3n+m} &=& \int_0^1 z^{3n+m-1} dz \\
\end{eqnarray*}
and summing the geometric series, we have the following triple integral
\begin{eqnarray*}
\int_0^1 \int_0^1 \int_0^1 \frac{z^3 dx dy dz}{(1-xz^3)(1-yz)}.
\end{eqnarray*}
Now doing the $x$ and $y$ integrations we have
\begin{eqnarray*}
I=\int_0^1  \frac{\ln(1-z) \ln(1-z^3)}{z} dz.
\end{eqnarray*}
Factorize the argument of the second logarithm ...
\begin{eqnarray*}
I= \underbrace{\int_0^1  \frac{\ln(1-z) \ln(1-z)}{z} dz}_{=2\zeta(3)} + \int_0^1  \frac{\ln(1-z) \ln(1+z+z^2)}{z} dz.
\end{eqnarray*}
So if you prefer my question is ... find a closed form for:
\begin{eqnarray*}
I_1 = - \int_0^1  \frac{\ln(1-z) \ln(1+z+z^2)}{z} dz.
\end{eqnarray*}
Integrating by parts gives:
\begin{eqnarray*}
I_1 = - \int_0^1  \frac{\ln(z) \ln(1+z+z^2)}{1-z} dz + \int_0^1  \frac{(1+2z)\ln(z) \ln(1-z)}{1+z+z^2} dz.
\end{eqnarray*}
and let us call these integrals $I_2$ and $I_3$ respectively.
All $3$ of these integrals are not easy for me to evaluate and any help with their resolution will be gratefully received.
 A: $$\boxed{I=\int_0^1  \frac{\ln(1-x) \ln(1-x^3)}{x}dx=\frac53\zeta(3) +\frac{2\pi^3}{27\sqrt 3} -\frac{\pi}{9\sqrt 3}\psi_1\left(\frac13\right)}$$
As mentioned in the question we have:
$$I=\int_0^1 \frac{\ln^2(1-x)}{x}dx+\int_0^1 \frac{\ln(1-x)\ln(1+x+x^2)}{x}dx=2\zeta(3)+J$$
We can make use of the following series:
$$ -\frac12 \ln(1-2x\cos t+x^2)=\sum_{n=1}^\infty \frac{\cos(nt)}{n} x^n,\quad |x|<1, t\in \mathbb R$$
$$\Rightarrow J=\int_0^1 \frac{\ln(1-x)\ln(1+x+x^2)}{x}dx=-2\sum_{n=1}^\infty \frac{\cos\left(\frac{2n \pi}{3}\right)}{n}\int_0^1 \ln(1-x) x^{n-1}dx $$
$$=2\sum_{n=1}^\infty \frac{\cos\left(\frac{2n \pi}{3}\right)}{n^2}H_n=2\Re \left(\sum_{n=1}^\infty \frac{z^n}{n^2}H_n\right),\quad z=e^{\frac{2\pi i}{3}}$$
Using the following generating function:
$$\sum_{n=1}^\infty \frac{x^n}{n^2}H_n=\operatorname{Li}_3(x)-\operatorname{Li}_3(1-x)+\operatorname{Li}_2(1-x)\ln(1-x)+\frac{1}{2}\ln x \ln^2(1-x)+\zeta(3)$$
And by plugging in the values found in this post yields the announced result, as we obtain:
$$J=\int_0^1 \frac{\ln(1-x)\ln(1-x+x^2)}{x}dx=\frac{2\pi^3}{27\sqrt 3}-\frac13\zeta(3) -\frac{\pi}{9\sqrt 3}\psi_1\left(\frac13\right)$$
