\begin{cases}
~\frac{\partial p}{\partial t}=c_0^2\rho_0\left(\frac{v_r}{r}+\frac{\partial v_r}{\partial r}\right)=c_0^2\rho_0\frac{1}{r}\frac{\partial}{\partial r}\left(rv_r\right)\\
~\frac{\partial v_r}{\partial t}=\frac{1}{\rho_0}\frac{\partial p}{\partial r}
\end{cases}
$$\frac{\partial^2 p}{\partial t^2}=c_0^2\rho_0\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial v_r}{\partial t}\right)=c_0^2\rho_0\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{1}{\rho_0}\frac{\partial p}{\partial r}
\right)$$
$$\frac{\partial^2 p}{\partial t^2}=c_0^2\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial p}{\partial r}
\right)=c_0^2\left(\frac{\partial p}{\partial r}
+\frac{1}{r}\frac{\partial^2 p}{\partial r^2}\right)$$
Search of particular solution with separated variables, on the form $p=R(r)T(t)$ :
$$\frac{T''}{T}=\frac{R'}{R}+\frac{1}{r}\frac{R''}{R}=\lambda=\text{constant}$$
$$\begin{cases}
T''-\lambda T=0\quad\implies\quad T(t)=C_1e^{-\sqrt{\lambda}\:t}
+C_2e^{\sqrt{\lambda}\:t}\\
R''+r R'-\lambda r R=0\quad\implies\quad\end{cases}$$
$ R(r)=e^{r(\lambda+\frac{r}{2})}\left(c_1 M\left( \frac{1-\lambda^2}{2}\:;\:\frac12\:;\:(\lambda+\frac{r}{2})^2\right) +c_2\:U\left( \frac{1-\lambda^2}{2}\:;\:\frac12\:;\:(\lambda+\frac{r}{2})^2\right) \right) $
$M$ denotes the Kummer confluent hypergeometric function of the first kind.
$U$ denotes the Kummer confluent hypergeometric function of the second kind.
http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html
http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheSecondKind.html
The general solution of the PDE is any linear combination of the particular solutions. On integral form :
$ P(r,t)=\int \left(F_1(\lambda)e^{-\sqrt{\lambda}\:t}
+F_2(\lambda)e^{\sqrt{\lambda}\:t} \right)e^{r(\lambda+\frac{r}{2})}\left(f_1(\lambda) M\left( \frac{1-\lambda^2}{2}\:;\:\frac12\:;\:(\lambda+\frac{r}{2})^2\right) +f_2(\lambda)\:U\left( \frac{1-\lambda^2}{2}\:;\:\frac12\:;\:(\lambda+\frac{r}{2})^2\right) \right) d\lambda$
$F_1$ , $F_2$ , $f_1$ , $f_2$ are arbitrary functions of $\lambda$. They have to be determined according to some initial and boundary conditions.
Alternatively the solution can be expressed on discret form instead of integral form, with a lot of arbitrary coefficients into it, again to be determined according to some initial and boundary conditions.
This is an arduous task requiring that the initial and boundary condition be very clearly and fully specified. Without a cautious definition of the conditions, they are an infinity of solutions for the PDE. With such confluent hyper geometric functions involved, one can fear very complicated calculus and solution. But simplifications are likely to occur for some particular cases, possibly avoiding the way of separation of variables and so avoiding the complicated special functions.
