I posted a similar question previously, here, but I realized I made a mistake in my original derivation of the equation I'm trying to determine a BC for so I'm going to try again. There were enough changes that I thought it'd better to make a new post instead of editing the old one a bunch.

I am looking to determine a boundary condition at $r=r_m$ for the following PDE:$$ \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)-\frac{2\mu_T}{r} \frac{\partial}{\partial r} \left(r \frac{\partial z(r,t)}{\partial r} \right) \frac{\partial z(r,t)}{\partial t}\tag{1} $$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 terms gives the following: $$ \left(\frac{\partial^2z}{\partial r^2}+\frac{1}{r} \frac{\partial z}{\partial r}\right)\left(\sigma_I+2\mu_T\frac{\partial z}{\partial t}\right)-\rho_Dgz=0 \tag{2}$$ $$\frac{\partial^2z}{\partial r^2}+\frac{1}{r} \frac{\partial z}{\partial r}-\frac{\rho_Dg}{\sigma_I+2\mu_T\frac{\partial z}{\partial t}}z=0 \tag{3}$$ I am trying to determine an analytical solution to either equation 2 or 3 so I can use it as a boundary condition so solve equation 1 numerically. Equation 3 is somewhat similar to the modified Bessel differential equation, but the $z$ coefficient throws that off by having $\frac{\partial z}{\partial t}$ in it.

In a simplified version of this problem $\mu_T=0$ thereby making the $z$ coefficient constant, and the second equation has 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)$ where $\lambda=\sqrt{\frac{\sigma_I}{\rho_Dg}}$. The modified Bessel function of the first kind is ignored as a possible solution since it is unbounded at lager argument values and $z(\infty,t)=0$ for this problem. 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$ $(r_m\approx 0.25\lambda)$ so this starts by using the asymptotic form of that Bessel function for smaller argument values, $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 matched to an analytical solution to the equation 1 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$, and the equation doesn't quite follow this previous solution method.

I've tried a few solutions but I'm unsure of their validity. First, I treat the nonlinear term like a constant so $\tau_v=-\frac{2\mu_T}{r} \frac{\partial}{\partial r} \left(r \frac{\partial z}{\partial r} \right) \frac{\partial z}{\partial t}$. With this, equation 1 then simplifies down to the modified Bessel function equal to the constant $\tau_v$. Going through the same process of matching the asymptotic form for smaller arguments to an analytical solution for smaller $r$ gives the boundary condition as $z(r_m,t)=\frac{1}{\sigma_I}\int_0^{r_m}{r(p-\tau_v)\mathrm{d}r}K_0\left(\frac{r_m}{\lambda} \right)-\frac{\tau_v(r_m)}{\rho_Dg}$.

The second method I've tried is using equation 3 but treating $\sqrt{\frac{\sigma_I+2\mu_T\frac{\partial z}{\partial t}}{\rho_Dg}}=\lambda_D$ as a constant to then say it has an analytical solution of $z=AK_0\left(\frac{r}{\lambda_D}\right)$. Once again, matching the asymptotic form to the analytical solution for small $r$ gives the boundary condition as $z(r_m,t)=\int_0^{r_m}{\frac{rp}{\sigma_I+2\mu_T\frac{\partial z}{\partial t}}\mathrm{d}r}K_0\left(\frac{r_m}{\lambda_D} \right)$. Since $\lambda_D$ is also a function of $r$ and $t$, I'm thinking to just use its value at $r=r_m$ when evaluating the Bessel function.

The third and final method I've tried is using separation of variables. So letting $z(r,t)=R(r)T(t)$ and substituting into the equation 1 with $p=0$: $$ \sigma_I\frac{\partial^2z}{\partial r^2}+\sigma_I\frac{1}{r} \frac{\partial z}{\partial r}=\rho_D gz-2\mu_T\frac{\partial z}{\partial t}\frac{\partial^2 z}{\partial r^2}-2\mu_T\frac{\partial z}{\partial t}\frac{1}{r} \frac{\partial z}{\partial r}\tag{4}$$ $$ \sigma_IT\frac{\partial^2R}{\partial r^2}+\sigma_I\frac{T}{r} \frac{\partial R}{\partial r}=\rho_D gTR-2\mu_TR\frac{\partial T}{\partial t}T\frac{\partial^2 R}{\partial r^2}-2\mu_TR\frac{\partial T}{\partial t}\frac{T}{r} \frac{\partial R}{\partial r}\tag{5}$$ 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(R\frac{\partial^2R}{\partial r^2}+\frac{R}{r}\frac{\partial R}{\partial r}\right)^{-1}=-\frac{2\mu_T}{\sigma_I} \frac{\partial T}{\partial t}=E\tag{6} $$ 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.

Once again, any advice or help you could provide would be greatly appreciated, thank you!