We can compute an asymptotic expansion of this sum using Mellin transforms. Let $f(x)$ be the base function $$f(x) = \frac{1}{3^x-1}.$$ The Mellin transform $f^*(s)$ of $f(x)$ is
$$ \mathfrak{M}(f(x); s) = f^*(s) = \int_0^\infty \frac{x^{s-1}}{3^x-1} dx
= \int_0^\infty \frac{1}{3^x} \frac{1}{1-3^{-x}} x^{s-1} dx \\ =
\int_0^\infty \frac{1}{3^x} \left(\sum_{q\ge 0} 3^{-qx} \right) x^{s-1} dx =
\sum_{q\ge 0} \int_0^\infty 3^{-(q+1)x} x^{s-1} dx \\=
\Gamma(s) \sum_{q\ge 0} \frac{1}{(\log 3)^s (q+1)^s} =
\frac{1}{(\log 3)^s} \Gamma(s) \zeta(s). $$
It follows that the Mellin transform of the harmonic sum
$$ g(x) = \sum_{n\ge 1} \frac{1}{3^{nx}-1}$$ is
$$ \mathfrak{M}(g(x); s) = g^*(s) = \frac{1}{(\log 3)^s} \Gamma(s) \zeta(s)^2.$$
Now invert to get the sum. We list the contributions from the main poles. Sum these to get the aymptotic expansion.
$$\begin{array}
\operatorname{Res}(g^*(x) x^{-s}; s=1) & = &
\frac{1}{\log 3} \left( (\gamma - \log \log 3) \frac{1}{x} - \frac{\log x}{x} \right)\\
\operatorname{Res}(g^*(x) x^{-s}; s=0) & = & \frac{1}{4} \\
\operatorname{Res}(g^*(x) x^{-s}; s=-1) & = & -{\frac {1}{144}}\,\log \left( 3 \right) x \\
\operatorname{Res}(g^*(x) x^{-s}; s=-3) & = &
-{\frac {1}{86400}}\,{x}^{3} \left( \log 3 \right) ^{3} \\
\operatorname{Res}(g^*(x) x^{-s}; s=-5) & = &
-{\frac {1}{7620480}}\,{x}^{5} \left( \log 3 \right)^5.
\end{array}$$
This partial expansion is already quite good, it gives $0.6821535092$ at $x=1$ whereas the precise value is $0.6821535026.$