Axial shortening of a torsion rod loaded with pure torque If you have a torsion rod and stress it with pure torsion you expect an axial shortening, too. If you have a lateral contraction property in your material (e.g. steel)
I anyhow cannot calculate the axial deformation caused by torsion e.g.
mu=0.3 % Poisson coefficient
sigma=0.1 %Radial strain in rad
l_0=10 %Inertial length of rod
E=210000 %Young's modulus

so what I want to calculate is
Delta_l=? %Axial displacement

According to St. Venant you won't have an axial displacement. But if you take the volume effects into account you will have...
 A: Let us see. The deformation in general can be seen as the superposition of bulk deformation (hydrostatic) plus shear deformation. In your case you only have the shear deformation which means volume is preserved.
So you need to find is the length of an outer fiber as it twists around along the final length of the rod $\ell$. The total twist is given by $\Theta = \frac{T \ell_0}{G\,J_T}$ which I assume you know how to calculate.
The helix has coordinates
$$\vec{r}(\theta) = \begin{pmatrix} R \cos\theta  & R \sin \theta & h \theta \end{pmatrix}$$
where $R$ is the radius and $h = \frac{\rm \ell}{\rm \Theta}$ the pitch of the helix. Note that the pitch decreases as the length shrinks. The length of the fiber is equal to the original length $\ell_0$ with:
$$ \ell_{0}=\int\,{\rm d}s=\int\,|{\rm d}\vec{r}|=\int_{0}^{\Theta}\sqrt{R^{2}+h^{2}}\;{\rm d}\theta\\=\Theta\sqrt{R^{2}+h^{2}}=\Theta\sqrt{R^{2}+\left(\frac{\ell}{\Theta}\right)^{2}}   $$
which yields the final length $\ell$ by
$$ \frac{ \ell_0}{ \Theta} = \sqrt{R^2 + \left(\frac{\ell}{\Theta}\right)^2} $$
$$ \boxed{ \ell = \sqrt{ \ell_0^2 - R^2 \Theta^2 } }$$
