I try to calculate $ \langle \boldsymbol{p} \rangle $ in the ground state $ \psi_{100} = \sqrt{\frac{Z^3}{\pi a^3}}\exp[-\frac{Z}{a}r] $
\begin{align} \langle \boldsymbol{p} \rangle = -i\hbar \frac{Z^{3}}{\pi a^{3}} \cdot 4\pi \int_{0}^{\infty} r^{2} \exp\left[-\frac{Z}{a}r\right] \left( \frac{\partial}{\partial r} \exp\left[-\frac{Z}{a}r\right] \right) \ {\rm d}r \ \boldsymbol{\mathrm{e}}_{r} = \frac{8Z}{a}i\hbar \boldsymbol{\mathrm{e}}_{r} \end{align}
I know it must be a wrong answer.
Yes, it is wrong. The expectation value of $\vec p$ is zero in any stationary state (with $H=\frac{{\vec p}^2}{2m}+V(\vec x)$). (Or any real-valued state whatsoever).
So, can we use the method in quantum to get the velocity $ v_n=\frac{Ze^2}{4\pi\epsilon_0\boldsymbol{\hbar}}\frac1n $?
If I were you, I would use the Virial theorem: $$ \langle 2T\rangle = m\langle \vec v^2\rangle = -\langle V\rangle = \frac{Z^2}{n^2}\;, $$ where I am using atomic units ($\hbar=m=e^2=a_0=1$).
The second equals sign follows from the Virial theorem, and the third equals sign follows from the definition of the potential $$ V = -\frac{Ze^2}{r}\;, $$ in Gaussian cgs units, and from the well-known expression $$ \frac{Z}{a_0 n^2} = \langle \frac{1}{r}\rangle $$ for a hydrogenic state.