I would like to solve a very simple case of 2D pressure wave propagation:
\begin{cases} ~\frac{\partial p}{\partial t}=c_0^2\rho_0\left(\frac{\partial v_x}{\partial x}+\frac{\partial v_y}{\partial y}\right)\\ ~\frac{\partial v_x}{\partial t}=\frac{1}{\rho_0}\frac{\partial p}{\partial x}\\ ~\frac{\partial v_y}{\partial t}=\frac{1}{\rho_0}\frac{\partial p}{\partial y} \end{cases}
Since this problem is actually radial, by using polar coordinates, I can rewrite the system as:
\begin{cases} ~\frac{\partial p}{\partial t}=c_0^2\rho_0\left(\frac{v_r}{r}+\frac{\partial v_r}{\partial r}\right)\\ ~\frac{\partial v_r}{\partial t}=\frac{1}{\rho_0}\frac{\partial p}{\partial r} \end{cases}
That I could have easily solved without the additional $v_r/r$ in the first equation. What is the method to consider?
\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.