I am trying to determine an outer boundary condition for the following PDE at $r=r_m$: $$ \frac{\sigma_I}{r} \frac{\partial}{\partial r} \left(r \frac{\partial z(r,t)}{\partial r} \right)=\rho_D gz(r,t)-p(r,t)-4 \mu_T \frac{\partial^2z(r,t)}{\partial r^2} \frac{\partial z(r,t)}{\partial t} $$ where $\sigma_I$, $\rho_D$, g, and $\mu_T$ are constants. The initial condition is $z(r,0)=0$ and the inner boundary condition at $r=0$ is that $\frac{\partial z}{\partial r}=0$. I know that at at $r\geq r_m$, $p=0$ for all $t$. Applying this and expanding the left side gives the following. $$ \frac{\partial^2z}{\partial r^2}+\frac{1}{r} \frac{\partial z}{\partial r}-\frac{z}{\lambda^2}=-\alpha\frac{\partial^2z}{\partial r^2} \frac{\partial z}{\partial t} $$ where $\lambda=\sqrt {\frac{\sigma_I}{\rho_Dg}}$ and $\alpha=\frac{4\mu_T}{\sigma_I}$. I am trying to determine an analytical solution to this equation to use as a boundary condition to solve the first equation numerically. The left hand side is the modified Bessel equation but obviously the right side changes that.

In a simplified version of this problem $\mu_T=0$, so $\alpha=0$ and the second equation becomes homogeneous with an analytical solution in the form of the modified Bessel function of the second kind of order zero, i.e. $z=AK_0\left(\frac{r}{\lambda} \right)$. I'm not sure why but the modified Bessel function of the first kind is ignored as a possible solution. The coefficient $A$ is determined by matching this solution for $r\geq r_m$ to the analytical solution for smaller $r$. $r_m$ is typically smaller than $\lambda$ so this starts by using the asymptotic form of that Bessel function, $z=AK_0\left(\frac{r}{\lambda} \right)\approx -A\left[ln\left(\frac{r}{2\lambda} \right)-\gamma_E\right]$, where $\gamma_E$ is the Euler constant. This asymptotic form is mathched to an analytical solution to the first equation at small $r$ where $\rho_Dgz$ is negligible but $p$ is no longer zero. Again with $\mu_T=0$ and those assumption, the first equation simplifies to $ \frac{\sigma_I}{r} \frac{\partial}{\partial r} \left(r \frac{\partial z}{\partial r} \right)=p$. Integrating once and applying the inner boundary condition gives $r\frac{dz}{dr}=-\frac{1}{\sigma_I}\int_0^{r_m}{rp\mathrm{d}r}$. The second integration gives $z=-\frac{1}{\sigma_I}\int_0^{r_m}{rp\mathrm{d}r}\ln(r/2\lambda)+...$ Matching coefficeints gives $A(t)=\frac{1}{\sigma_I}\int_0^{r_m}{rp\mathrm{d}r}$. Thus the outer boundary condition used is $z(r_m,t)=\frac{1}{\sigma_I}\int_0^{r_m}{rp\mathrm{d}r}K_0\left(\frac{r_m}{\lambda} \right)$. I'm hoping to do some similar analysis with the second equation shown here where $\mu_T\neq0$, thus $\alpha\neq0$ and the equation doesn't quite follow this previous solution method.

I have tried a few things so far but with no success. First, I pretend the right side is some constant $p_v=-\alpha\frac{\partial^2z}{\partial r^2} \frac{\partial z}{\partial t}$ and follow the same method I just described for the simplified case. This gives the boundary condition as $z(r_m,t)=\frac{1}{\sigma_I}\int_0^{r_m}{r(p-p_v)\mathrm{d}r}K_0\left(\frac{r_m}{\lambda} \right)-\frac{p_v(r_m)}{\rho_Dg}$. Unfortunately, the numerical solution fails to meet any integration tolerances with this.

Secondly, I've tried matching the second equation shown here to an alternate form of the Bessel equation which was given in a 1958 book by Bowman (and also on Wolframalpha). The alternate form looks like this: $$ \frac{\mathrm{d}^2z}{\mathrm{d}r^2}+\frac{1-2a}{r} \frac{\mathrm{d}z}{\mathrm{d}r}+\left(b^2c^2r^{2c-2}+\frac{a^2-n^2c^2}{r^2}\right)z=0 $$ and has the following solutions $$ z= \begin{cases} r^a\left[AJ_n(br^c)+BY_n(br^c)\right] &\text{ for integer }n \\ \\ r^a\left[AJ_n(br^c)+BJ_{-n}(br^c)\right] &\text{ for noninteger }n \end{cases} $$ Since the $z$ coefficient in my second equation has no $r$ dependence, $2c-2=0$ so $c=1$, and $a^2-n^2c^2=0$ so $a=n$. Thus, the alternate form simplifies down to: $$ \frac{\mathrm{d}^2z}{\mathrm{d}r^2}+\frac{1-2a}{r} \frac{\mathrm{d}z}{\mathrm{d}r}+b^2z=0 $$ and the simplified solutions are: $$ z= \begin{cases} r^a\left[AJ_a(br)+BY_a(br)\right]&\text{ for integer }n \\ \\ r^a\left[AJ_a(br)+BJ_{-a}(br)\right]&\text{ for noninteger }n \end{cases} $$ Rewriting my second equation into this alternate simplified form gives: $$ \frac{\partial^2z}{\partial r^2}+\frac{1}{1+\alpha \frac{\partial z}{\partial t}}\frac{1}{r} \frac{\partial z}{\partial r}-\frac{1}{\lambda^2\left(1+\alpha \frac{\partial z}{\partial t}\right)}z=0 $$ Treating the coefficients as constants, I write the general solution as: $$ z= \begin{cases} r^a\left[AJ_a(-ibr)+BY_a(-ibr)\right] & \text{ for integer }n\\ \\ r^a\left[AJ_a(-ibr)+BJ_{-a}(-ibr)\right] & \text{ for noninteger }n \end{cases} $$ where $1-2a=\frac{1}{1+\alpha \frac{\partial z}{\partial t}}$ and $b^2=\frac{1}{\lambda^2\left(1+\alpha \frac{\partial z}{\partial t}\right)}$. From here, I'm not really sure where to go or if the analysis up to this point is even correct. Since the arguments of the Bessel functions have the imaginary unit, $i$, maybe I can switch over to the modified Bessel functions. Still need a way to get the $A$ and $B$ however. Possibly ignore one of the Bessel functions like was done in the previous solution method and then match the coefficient to an analytical solution for smaller $r$ where $\rho_Dgz$ is negligible but $p$ is not.

The third and final method I've tried is using separation of variables. So letting $z(r,t)=R(r)T(t)$ and substituting in to the second equation shown here gives: $$ T\frac{\partial^2R}{\partial r^2}+\frac{T}{r} \frac{\partial R}{\partial r}-T\frac{R}{\lambda^2}=-\alpha T\frac{\partial^2R}{\partial r^2}R\frac{\partial T}{\partial t} $$ The $T$ divides out and this can then be separated to get the following: $$ \left(\frac{\partial^2R}{\partial r^2}+\frac{1}{r} \frac{\partial R}{\partial r}-\frac{R}{\lambda^2}\right)\left(\frac{\partial^2R}{\partial r^2}R\right)^{-1}=-\alpha \frac{\partial T}{\partial t}=E $$ where $E$ is some constant. Unfortunately this still leaves two issues: (1) it doesn't physically make sense for the temporal derivative to be equal to a constant. $z$ is describing a surface which moves under an applied pressure. Imagine bouncing a basketball at the center of a trampoline and $z$ is describing the axisymmetric surface shape of the trampoline. So every time the basketball hits the trampoline, the surface deforms slightly and then returns to its original flat shape as the ball rebounds away from the surface. The surface velocity in the $z$ direction can't be constant or it's like your trampoline surface is constantly moving in one direction. Ignoring the physics for a second, issue (2) is the spatial ODE is still a second order nonlinear equation which I'm unsure of how to solve analytically. If that could be done, I could possibly match this analytical solution for $r=r_m$ to an analytical solution for smaller $r$ to determine the value of the constant $E$. But of course that doesn't reconcile the first issue.

Any advice or help you could provide would be greatly appreciated, thank you!