I am currently struggling with the same problem, getting a bit further but not all the way. If you solved it in the meantime, I'd be glad to hear about it.
EDIT 2: Improved answer some more as my technique improved some more.
EDIT 3: The pull request (linked below) now has somewhat better implementation then discussed here that reduces the error a bit further, but still breaks for some values.
C++ code (written for Stan math library) can be found at https://github.com/stan-dev/math/pull/1121
Case 0: $x$ small relative to $\alpha$
I use this when $\alpha > 15$ AND $10x < \sqrt{\alpha + 1}$. The formula can be found on Wikipedia on Bessel K and also as formula 1.10 of Temme, Journal of Computational Physics, vol 19, 324 (1975). It has:
$$
K_\alpha(x) \simeq \frac{\Gamma(\alpha)}{2} \left( \frac{2}{x} \right)^\alpha
$$
Case 1: Small to medium $\alpha$ and $x$
My main workhorse is using the logarithm of Equation 26 of Rothwell: Computation of the
logarithm of Bessel functions of complex argument and fractional order.
which has (maintaining their notation)
$$
K_\nu(z) = \frac{\sqrt{\pi}}{\Gamma(\nu + \frac{1}{2})}(2z)^{-\nu}e^{-z}
\int_0^1 [\beta e^{-u^\beta}(2z + u^\beta)^{\nu-1/2} u^{n-1} + e^{-1/u} u^{-2\nu-1}(2zu + 1)^{\nu-1/2}] \mathrm{d}u
$$
Where $\beta = \frac{2n}{2\nu+1}$ and the authors suggest using $n = 8$.
You still can't push the log through the integral, but for large areas of the parameter space, it is OK to just compute the integral and take the log afterwards.
The integral ceases to be numerically tractable around roughly $\nu > 50$ or when $\nu > 0.5$ AND $log(z) > \frac{300}{\nu - 0.5} - log(2)$ (one of the summands turn to infinity).
Fruther, the formula ceases to be very accurate for small $z$ and small $\nu$, where below roughly $z < 10^-4$ the relative error starts climbing from around $10^{-30}$ to around $10{^-2}$ for $z < 10^{-7}$
There is a C++ implementation of the Rothwell formula at
https://github.com/stan-dev/stan/wiki/Stan-Development-Meeting-Agenda/0ca4e1be9f7fc800658bfbd97331e800a4f50011
Case 2: $\alpha \gg x$
Here the log of the approximate formula you've written works OK, but I got a bit better results with the same formula as in Case 0.
Case 3: $x \gg \alpha$
Asymptotic formula such as 10.4.2 here: https://dlmf.nist.gov/10.40 can be used. Once again, you need to compute the sum on non-log scale (it has negative elements), and take the log afterward, but this has worked great for me. Summing up to 10 first elements worked OK in my case.
Remaining cases
For large $x$ with comparably large $\alpha$, I still haven't found a reliable solution.
Here is a plot of error with respect to recursive formula for consecutive neighbours with the approach I've described (ratio is the relative error of the actual value, I also test the gradient of this function wrt. both parameters as computed by Stan's autodiff - the gradient is a tougher problem, but probably not relevant here):
The biggest relative error shown is 3.8e-02 for $\alpha = 148$ and $x = 105$.
Best Answer
By $(10.22.49)$ and $(15.4.19)$ in the DLMF, we find \begin{align*} & \int_0^{ + \infty } {\frac{{t^3 }}{{\sinh t}}J_1 (at)dt} = 2\int_0^{ + \infty } {t^3 e^{ - t} \frac{1}{{1 - e^{ - 2t} }}J_1 (at)dt} \\ & = 2\int_0^{ + \infty } {t^3 \sum\limits_{n = 0}^\infty {e^{ - (2n + 1)t} } J_1 (at)dt} \\ & = 2\sum\limits_{n = 0}^\infty {\int_0^{ + \infty } {t^3 } e^{ - (2n + 1)t} J_1 (at)dt} \\ & = 24a\sum\limits_{n = 0}^\infty {\frac{1}{{(2n + 1)^5 }}F\!\left( {\frac{5}{2},3;2; - \frac{{a^2 }}{{(2n + 1)^2 }}} \right)} \\ & = 24a\sum\limits_{n = 0}^\infty {\left( {(2n + 1)^2 - \frac{{a^2 }}{4}} \right)\frac{1}{{((2n + 1)^2 + a^2 )^{7/2} }}} \\ &= 24a\sum\limits_{n = 0}^\infty {\frac{1}{{((2n + 1)^2 + a^2 )^{5/2} }}} - 30a^3 \sum\limits_{n = 0}^\infty {\frac{1}{{((2n + 1)^2 + a^2 )^{7/2} }}} \end{align*} provided $\Re a>0$. Here $F$ stands for the hypergeometric function.
Addendum. I shall give the asymptotics for $a\to +\infty$. Let us introduce the generalised Mathieu series via $$ S_{\mu ,\gamma } (a;\lambda ) = \sum\limits_{n = 1}^\infty {\frac{{n^\gamma }}{{(n^\lambda + a^\lambda )^\mu }}} \quad \quad (\mu > 0,\quad \lambda > 0,\quad \lambda \mu - \gamma > 1). $$ With this notation $$ \sum\limits_{n = 0}^\infty {\frac{1}{{((2n + 1)^2 + a^2 )^\mu }}} = S_{\mu ,0} (a;2) - 4^{ - \mu } S_{\mu ,0} (a/2;2) $$ provided $\mu>\frac{1}{2}$. The precise asymptotics of the generalised Mathieu series was derived in this paper. In particular, $$ S_{\mu ,0} (a;2) = \frac{{\sqrt \pi \Gamma \left( {\mu - \frac{1}{2}} \right)}}{{2\Gamma (\mu )a^{2\mu - 1} }} - \frac{1}{{2a^{2\mu } }} + \frac{{\pi ^\mu }}{{\Gamma (\mu )a^\mu }}e^{ - 2\pi a} (1 + o(1)) $$ as $a\to +\infty$. Consequently, $$ \sum\limits_{n = 0}^\infty {\frac{1}{{((2n + 1)^2 + a^2 )^\mu }}} = \frac{{\sqrt \pi \Gamma \left( {\mu - \frac{1}{2}} \right)}}{{4\Gamma (\mu )a^{2\mu - 1} }} - \frac{{\pi ^\mu }}{{\Gamma (\mu )(2a)^\mu }}e^{ - \pi a} (1 + o(1)) $$ as $a\to +\infty$. Combining this with the exact series representation above, we find $$ \int_0^{ + \infty } {\frac{{t^3 }}{{\sinh t}}J_1 (at)dt} = \frac{{\sqrt 2 \,\pi ^3 }}{{a^{1/2} }}e^{ - \pi a} (1 + o(1)) $$ as $a\to +\infty$. More precise asymptotics can be derived by using more terms in the asymptotic expansion of the generalised Mathieu series which can be found in the paper cited above.