How to combine thermal velocity concept and special relativity?

The assumption that the thermal energy is $\frac{3}{2}k_bT$ is actually only valid at non-relativistic temperatures. In general we have to use the equipartition theorem to find the relation between temperature and energy:

\begin{equation} \left< x_m \frac{\partial E}{\partial x_n} \right> = \delta_{mn} k_BT, \end{equation} where $x$ can be a coordinate or conjugate momentum.

Just taking the one-dimensional case for simplicity, in the Newtonian regime, $E = \frac{mv^2}{2}$, so that $v=\sqrt{\frac{k_BT}{m}}$. But in the relativistic case, $E = \sqrt{p^2c^2 + m_0^2c^4}$. This means that

\begin{equation} \frac{c^2p^2}{\sqrt{p^2c^2 + m_0^2c^4}} = k_BT, \end{equation} so \begin{equation} p^2 = \frac{k_B^2T^2c^2 \pm \sqrt{k_B^4T^4c^4 + 4k_B^2T^2m_0^2c^8}}{2c^4}. \end{equation}

As $T \rightarrow \infty$ we get $p = k_BT/c$, or $E=k_BT$ (as the mass-term in the energy becomes negligible compared to the momentum). This is the well-known energy-equation for the ultra-relativistic gas. As $T \rightarrow 0$ we get $v = \sqrt{\frac{k_BT}{m_0}}$, the Newtonian result.

What you are looking for is the mean of the Maxwell–Jüttner distribution, $$f(\gamma)=\frac{\gamma^2\beta}{\theta K_2(1/\theta)}e^{-\gamma/\theta}$$ where $\beta=v/c=\sqrt{1-1/\gamma^2}$, $\theta=k_BT/mc^2$ and $K_2$ is the Bessel function of the second kind.

So the expected $\gamma$ would be $$E[\gamma]=\frac{1}{\theta K_2(1/\theta)}\int_1^\infty \gamma^3 \left (\sqrt{1-1/\gamma^2} \right )e^{-\gamma/\theta}d\gamma$$ Unfortunately there doesn't seem to be any closed form formula for it. Here is a plot of my numeric results in Matlab:

Here is a derivation of anisotropic Maxwell–Jüttner distributions, which may or may not be useful.