Determine PDE Boundary Condition via Analytical solution

[tse8682]

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}=-\frac{4\mu_T}{\sigma_I}\frac{\partial^2z}{\partial r^2} \frac{\partial z}{\partial t} $$ where $\lambda=\sqrt {\frac{\sigma_I}{\rho_Dg}}$. 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$, 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)$. The modified Bessel function of the first kind is ignored as a possible solution because it becomes unbounded at larger argument values. 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$, $\left(r_m\approx 0.25\lambda\right)$, 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 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 is done via integration by parts and gives $z=-\frac{1}{\sigma_I}\int_0^{r_m}{rp\mathrm{d}r}\ln(r/2\lambda)+...$ some other stuff. 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 have tried a few things so far but with no success. First, I pretend the nonlinear term is some constant $\varepsilon_v=-4\mu_T\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-\varepsilon_v)\mathrm{d}r}K_0\left(\frac{r_m}{\lambda} \right)-\frac{\varepsilon_v(r_m,t)}{\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 any }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 any }a \\
\\
r^a\left[AJ_a(br)+BJ_{-a}(br)\right]&\text{ for noninteger }a
\end{cases} $$ Rewriting my second equation into this alternate simplified form gives:
$$ \frac{\partial^2z}{\partial r^2}+\frac{1}{1+\frac{4\mu_T}{\sigma_I}\frac{\partial z}{\partial t}}\frac{1}{r} \frac{\partial z}{\partial r}-\frac{1}{\lambda^2\left(1+\frac{4\mu_T}{\sigma_I}\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(-irb)+BY_a(-irb)\right] & \text{ for any }a\\
\\
r^a\left[AJ_a(-irb)+BJ_{-a}(-irb)\right] & \text{ for noninteger }a
\end{cases}
$$ or $$
z=
\begin{cases}
r^a\left[AI_a(br)+BK_a(br)\right] & \text{ for any }a\\
\\
r^a\left[AI_a(br)+BI_{-a}(br)\right] & \text{ for noninteger }a
\end{cases}
$$
where $1-2a=\frac{1}{1+\frac{4\mu_T}{\sigma_I}\frac{\partial z}{\partial t}}$ and $b^2=\frac{1}{\lambda^2\left(1+\frac{4\mu_T}{\sigma_I}\frac{\partial z}{\partial t}\right)}$. Again, since $I_a$ is unbounded at larger arguments, $A$ can be asumed to be zero. Technically, $a$ and $b$ are not constants but functions of both $r$ and $t$. I'm hoping I can just treat them as temporal functions since I'm only concerned with a single location, $r=r_m$. I'm unsure if this is a reasonable assumtion though. From here, I've been trying to determine an analytical solution for smaller values of $r$ to determine the value of $B$ as was done in the simplified verion of this problem. At smaller $r$, $\rho_Dgz$ is negligible but $p$ is not so the first equation simplifies to:
$$ \frac{\sigma_I}{r} \frac{\partial}{\partial r} \left(r \frac{\partial z}{\partial r} \right)=-p-4\mu_T\frac{\partial^2z}{\partial r^2} \frac{\partial z}{\partial t} $$ From here, I've looked at 2 options to find the inner analytical solution: (1) treat the nonlinear term as a constant again, $\varepsilon_v=-4\mu_T\frac{\partial^2z}{\partial r^2} \frac{\partial z}{\partial t}$ which gives me $B=\frac{1}{\sigma_I}\int_0^{r_m}{r(p-\varepsilon_v)\mathrm{d}r}$. Unfortunately the numerical solution does not converge using this method either. Option (2) is to combine the nonlinear term with the left side. This gets a little messy but looks like this. Expanding the left side and pulling out some variables:
$$ \frac{\partial^2z}{\partial r^2}+\frac{1}{1+\frac{4\mu_T}{\sigma_I}\frac{\partial z}{\partial t}}\frac{1}{r} \frac{\partial z}{\partial r}=-\frac{p}{\sigma_I\left(1+\frac{4\mu_T}{\sigma_I}\frac{\partial z}{\partial t}\right)} $$ $$ \frac{\partial^2z}{\partial r^2}+(1-2a)\frac{1}{r} \frac{\partial z}{\partial r}=-\frac{p}{\sigma_I}(1-2a) $$ $$ \frac{1}{r^{1-2a}} \frac{\partial}{\partial r} \left(r^{1-2a} \frac{\partial z}{\partial r} \right)=-\frac{p}{\sigma_I}(1-2a)$$ $$ \frac{\partial z}{\partial r}=-\frac{r^{-1+2a}}{\sigma_I}\int_0^r{p(1-2a)r^{1-2a}dr}$$ $$
z=
\begin{cases}
-\frac{1}{\sigma_I}\int_0^r{prdr}\ln{r}+...& \text{ for }a=0\\
\\
-\frac{r^{2a}}{\sigma_I2a}\int_0^r{p(1-2a)r^{1-2a}dr}+\frac{1}{\sigma_I}\int_0^r{p(1-2a)\frac{r}{2a}dr}+C_1 & \text{ for }a\neq0
\end{cases}
$$ I'm hoping to match this solution to the asymptotic form the bessel equation for small $r$ somehow. Those forms look like this:
$$
BK_a\approx
\begin{cases}
-B\ln{\frac{r}{2}}-B\gamma_E& \text{ for }a=0\\
\\
B\frac{\Gamma(a)}{2}\left(\frac{2}{r}\right)^a & \text{ for }a\gt0
\end{cases}
$$ The first case is simple to match the coefficients and get $B$ as its the same as the simplified version, but I have not been able to figure out a way to match the coefficeint for the case that $a\neq0$. This is a rather complex possible solution and would be difficult to implement even if I did find a way to determine $B$ for the other case.

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 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}=-\frac{4\mu_T}{\sigma_I} 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}=-\frac{4\mu_T}{\sigma_I}\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 the analytical solution for smaller $r$ to determine the value of the constant $E$. But of course, that doesn't reconcile the first issue.

If you've made it this far, thanks for just reading through the whole thing. Any advice or help you could provide would be greatly appreciated, thank you!

 

[tse8682]

I realized I made a mistake in the derivation of the original equation. Instead it should read: $$ \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} $$ Applying the $p=0$ boundary condition and rearranging then gives:$$ \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}$$ This is even closer to the modified Bessel function than I originally thought, except the $z$ coefficient is dependent on $r$ and $t$ as opposed to being a constant since $\frac{\partial z}{\partial t}$ is there.

I've tried a few solutions to this new equation as well but I'm unsure of their validity. First, I again 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}$ where $\lambda=\sqrt{\frac{\sigma_I}{\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 again 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 and $\lambda=\sqrt{\frac{\sigma_I}{\rho_Dg}}$. Unfortunately this still leaves the same two issues I described in the original post: (1) it doesn't make sense physically and (2) still leaves a rather complicated second order nonlinear ODE to solve.

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