I wish to evaluate $\sum\limits_{n\geqslant1}\frac1n\sum\limits_{d\mid n}\frac{d}{n^2+d}.$
Some observations:
Let $f(n)=\sum\limits_{d\mid n}\frac{d}{n^2+d}$.
Then $f(p)=\frac{p^2+p+2}{(p^2+1)(p+1)},\ f(2)=\frac8{15},\ f(4)=\frac{283}{765}.\ \ f$ isn't multiplicative, there is no way to relate $f(4)$ to $f(2),$ and trying to add the summation for $f$ to itself by considering $\frac{d}{n^2+d}+ \frac{n\mid d}{n^2+n\mid d}= \frac{n^2+d^2n+2d}{(n^2+d)(nd+1)}$ doesn't look promising.
Because none of these approaches worked, I turned to changing the order of summation.
Each summand for $f(n)$ is $\leqslant \frac{n}{n^2+n}=\frac1{n+1}<\frac1n$ so, the sum is $\leqslant\sum\limits_n \frac{d(n)}{n^2}\leqslant\sum\limits_n \frac{2\sqrt{n}}{n^2},$ which converges. Each term is positive, so the sum absolutely converges and we may rearrange the order of summation to sum over $d$ and $n=dk$ for $k\geqslant 1.$ We get $\sum\limits_{d\geqslant1} \sum\limits_{k\geqslant1}\frac1{dk} \frac{d}{(dk)^2+d}=\sum\limits_{d\geqslant 1}\frac1{d}\sum\limits_{k\geqslant 1}\frac1{k(dk^2+1)}=\sum\limits_{d \geqslant1}\frac1d\sum\limits_{k\geqslant1}\left(\frac1k- \frac{dk}{dk^2+1}\right).$ But now we have a difference of $2$ divergent series. How do we proceed when the partial fraction decomposition turns against you? Even the sum for $d=1$ does not appear to have a nice closed form.
Update: The problem has a typo and should have $d^2+n$ in the denominator. I found a solution to the corrected version.
Let $f(n) = \sum\limits_{d|n} \frac{d}{n+d^2}$ so that we seek $\sum\limits_{n \ge 1} \frac{f(n)}{n}.$ Since $\frac{n/d}{n+(n/d)^2} = \frac{d}{n+d^2},$ we have $f(n) \le 2\sum\limits_{d|n \atop d \ge \sqrt{n}} \frac{d}{n+d^2} < 2\sum\limits_{d|n \atop d \ge \sqrt{n}} \frac{1}{d} = \frac{2}{n} \sum\limits_{d|n \atop d \le \sqrt{n}} d \le \frac{d(n)\sqrt{n}}{n} < \frac{1}{n^{1/4}}$ for $n$ sufficiently large since $d(n) = o(n^{\epsilon})$ for any $\epsilon>0$ (here we chose $\epsilon = 1/4$)*. By comparison to $\sum\limits_{n \ge 1} n^{-5/4},$ the sum converges. Each term is positive, so the sum absolutely converges and we may rearrange the order of summation to sum over $d$ and $n = dk$ for $k \ge 1.$ We get $\sum\limits_{d \ge 1} \sum\limits_{k \ge 1} \frac{1}{dk} \frac{d}{dk+d^2} = \sum\limits_{d \ge 1} \frac{1}{d} \sum\limits_{k \ge 1} \frac{1}{k(d+k)} = \sum\limits_{d \ge 1} \frac{1}{d^2} \sum\limits_{k \ge 1} \left(\frac{1}{k} - \frac{1}{d+k}\right) = \sum\limits_{d \ge 1} \frac{H_d}{d^2}.$
Using the integral representation $H_n = \int_0^1 \frac{1-x^n}{1-x} \, dx,$ we get that the sum is $\int_0^1 \frac{g(x)}{1-x} \, dx$ where $g(x) = \sum\limits_{n \ge 1} \frac{1-x^n}{n^2}.$ Since $g(1) = 0$ and $xg'(x) = -\sum\limits_{n \ge 1} \frac{x^n}{n} = \log(1-x),$ we get $g(x) = -\int_x^1 \frac{\log(1-t)}{t} \, dt.$ Using integration by parts we get $\int_0^1 \frac{g(x)}{1-x} \, dx = -\log(1-x)g'(x) \big|_0^1 + \int_0^1 g'(x)\log(1-x) \, dx = \int_0^1 \frac{\log(1-x)^2}{x} \, dx = \int_0^1 \frac{\log(x)^2}{1-x} \, dx.$ Let $I_n = \int_0^1 \log(x)^2 x^n \, dx.$ By expanding $\frac{1}{1-x}$ as a power series and exchanging the order of summation and integration, we get $\int_0^1 \frac{\log(x)^2}{1-x} = \sum\limits_{n \ge 0} I_n.$
We turn our attention to $I_n = \int_0^1 \log(x)^2 x^n \, dx = \int_{-\infty}^0 u^2 e^{u(n+1)} \, du = \frac{1}{(n+1)^3} \int_{-\infty}^0 v^2 e^v \, dv = \frac{1}{(n+1)^3} \int_0^{\infty} v^2 e^{-v} = \frac{2!}{(n+1)^3}.$ Finally, the original sum is nothing but $2 \sum\limits_{n \ge 0} \frac{1}{(n+1)^3} = 2\zeta(3).$
*The elementary bound $d(n) \le 2\sqrt{n}$ is almost but not quite enough (you get the harmonic series), so we must resort to a result of Apostol. If there is a very fast proof of $d(n) = o(n^c)$ for some $c<1/2,$ I'd like to hear it.
By switching the order of summation and by applying the inverse Laplace transform we have $$S=\sum_{n\geq 1}\sum_{d\mid n}\frac{d}{n(n^2+d)} = \sum_{d\geq 1}\sum_{k\geq 1}\frac{1}{k(k^2 d^2+d)}=\int_{0}^{+\infty}\frac{ds}{e^s-1}\sum_{k\geq 1}\frac{1-e^{-s/k^2}}{k}\,ds=\sum_{k\geq 1}\frac{H_{1/k^2}}{k} $$ and now we may invoke the Maclaurin series of $H_s$ in a right neighbourhood of the origin $$ H_s = \zeta(2)s-\zeta(3)s^2+\zeta(4)s^3-\ldots $$ in order to get $$ S = \frac{1}{2}+\sum_{m\geq 2}(-1)^m\left[\zeta(m)\zeta(2m-1)-1\right] $$ which turns $S$ into an essentially geometric series, whose value is $\approx 1.29534$.