Podemos calcular una expansión asintótica de esta suma utilizando las transformadas de Mellin. Sea $f(x)$ sea la función base $$f(x) = \frac{1}{3^x-1}.$$ La transformada Mellin $f^*(s)$ de $f(x)$ es $$ \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). $$ Se deduce que la transformada de Mellin de la suma armónica $$ g(x) = \sum_{n\ge 1} \frac{1}{3^{nx}-1}$$ es $$ \mathfrak{M}(g(x); s) = g^*(s) = \frac{1}{(\log 3)^s} \Gamma(s) \zeta(s)^2.$$ Ahora invierte para obtener la suma. Enumeramos las contribuciones de los polos principales. Suma de estos para obtener la expansión asintótica. $$\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}$$
Esta ampliación parcial ya es bastante buena, da $0.6821535092$ en $x=1$ mientras que el valor exacto es $0.6821535026.$