Does the following hold?: $$ \sum_{a, b \in \mathbb{N}^+, \ \gcd(a,b) = 1} \frac{1}{ab(a+b)} \ = \ 2 $$ Numerical computations suggest this may hold, but on the other hand it would be quite surprising if that sum is indeed an integer.
If it indeed holds: is this a known identity?
This is a known identity: it is stated explicitly under Corollary 1 in Tornheim: Harmonic double series, Amer. J. Math. 72 (1950), 303-314. See also (10) in Mordell: On the evaluation of some multiple series, J. London Math. Soc. 33 (1958), 368-371, and see this MSE post for some neat derivations. In fact many other references are available, just search for "Tornheim sum" in Google and MathSciNet.
Added. The other responses, which were made before this addition, make it clear that the statement in the original post can be reduced to $\zeta(2,1)=\zeta(3)$. This famous identity on multiple zeta values was discovered by Euler (1775); see this MO post. A generalization $$\sum_{\substack{k_1+\dotsb+k_r=k\\k_1\geq 2;\ k_2,\dotsc,k_r\geq 1}}\zeta(k_1,k_2,\dots,k_r)=\zeta(k),\qquad k\geq r+1,$$ was conjectured by Moen (1990) and proved by Granville (1996). For a transparent proof of this general sum formula, see Seki-Yamamoto (2019) or its arXiv version.
I have no idea if it is known, but here is the proof:
First of all, we can remove the coprimality condition by factoring out the gcd and just prove that $S=\sum_{a, b = 1}^\infty \frac{1}{ab(a+b)} = 2\zeta(3)$. We will first fix $a$ and take the sum for $b$ from $1$ to infinity. $$\sum_{b = 1}^\infty \frac{1}{ab(a+b)} = \frac{1}{a^2}\sum_{b = 1}^\infty \left(\frac{1}{b} - \frac{1}{b+a}\right) = \frac{1}{a^2}\left(1 + \frac{1}{2} + \ldots + \frac{1}{a}\right).$$
So, we get $$S = \sum_{a = 1}^\infty \frac{1}{a^2}\sum_{k = 1}^a \frac{1}{k}.$$
We will change the variables to $a = k + s$ so now $k$ goes from $1$ to infinity and $s$ goes from $0$ to infinity. Term with $s = 0$ will give us $\zeta(3)$. The rest is $$S_1 = \sum_{k = 1}^\infty \sum_{s = 1}^\infty \frac{1}{k(k+s)^2}$$ and $S = S_1 + \zeta(3)$. Obviously, $$S_1 = \sum_{k = 1}^\infty \sum_{s = 1}^\infty \frac{1}{s(k+s)^2}$$ (we just switched the roles of $k$ and $s$). Adding these two expressions for $S_1$ we get $$2S_1 = \sum_{k = 1}^\infty \sum_{s = 1}^\infty \frac{1}{(k+s)^2}\left(\frac{1}{k}+\frac{1}{s}\right) = \sum_{k = 1}^\infty \sum_{s = 1}^\infty \frac{1}{ks(k+s)} = S.$$
So, $2S_1 = S$ and $S = S_1 + \zeta(3)$, hence $S = 2\zeta(3)$.
I'll use Iverson notation. Since $\mathrm{gcd}(a,b)=\mathrm{gcd}(a,a+b)$, the sum is $$S=\sum_{(a,b)\in S}\left(\frac1a+\frac1b\right)\frac1{(a+b)^2}=2\sum_{(a,b)\in S}\frac1a\cdot\frac1{(a+b)^2}=2\sum_{0<a<c}\frac{\varepsilon_{a,c}}{ac^2},$$ where $\varepsilon_{a,c}=[\mathrm{gcd}(a,c)=1]$. Therefore, $\varepsilon_{a,c}=\sum_{d\mid\mathrm{gcd}(a,c)}\mu(d)=\sum_{d\in\mathbb N^+}[d\mid a][d\mid c]\mu(d)$ and thus, $$\frac S2=\sum_{d\in\mathbb N^+}\sum_{0<a<c}\frac{[d\mid a]}a\frac{[d\mid c]}{c^2}\mu(d)=\sum_{d\in\mathbb N^+}\sum_{0<a_0<c_0}\frac{\mu(d)}{(da_0)(dc_0)^2}=\sum_{d\in\mathbb N^+}\frac{\mu(d)}{d^3}\sum_{0<a_0<c_0}\frac1{a_0c_0^2}=\frac{\zeta(1,2)}{\zeta(3)}=1.$$
