# Heat equation in cylindrical coordinates with Neumann boundary condition

Given a cylinder of internal radius $r_0$ and external radius $r_1$, the heat equation in cylindrical coordinates that represents the behaviour of the temperature inside the cylinder, can be written as: $$\frac{dT(\alpha,r,t)}{dt}=D[\frac{d^2 T(\alpha,r,t)}{dr^2}+\frac{1}{r}\frac{dT(\alpha,r,t)}{dr}+\frac{1}{r^2}\frac{d^2T(\alpha,r,t)}{d\alpha^2}+\frac{d^2T(\alpha,r,t)}{dz^2}]\qquad1$$ If we suppose the temperature isn't depending by $\alpha$ and $z$, the previous equation can be simplified to: $$\frac{dT(r,t)}{dt}=D[\frac{d^2 T(r,t)}{dr^2}+\frac{1}{r}\frac{dT(r,t)}{dr}]\qquad2$$ Assuming the following boundary and initial conditions: $$T(r,0)=T_0$$ for $r_0\le r \le r_1$ $$\partial_r T(r,t)=T_1sin(\omega t)$$ calculated in $r_1$, how can be found the solution of the equation $(2)$?

-
You have $T_{rr}$, a second derivative in space, so this problem needs two boundary conditions. – in_wolframAlpha_we_trust Oct 10 '12 at 13:56
@in_wolfram_we_trust: I see, but unfortunately, I have only this. – Riccardo.Alestra Oct 10 '12 at 14:35
$\partial_r T(r,t)=T_1\sin\omega t$ ? That means you can just find $T(r,t)$ by integration with respect to $r$ ? Is this a typo? It should be modify to $\partial_r T(r_0,t)=\partial_r T(r_1,t)=T_1\sin\omega t$ ? – doraemonpaul Oct 11 '12 at 5:16
@doraemonpaul: In fact I wrote $\partial_r T(r,t)$ calculated in $r_1$ – Riccardo.Alestra Oct 11 '12 at 8:41
@Riccardo.Alestra: I don't understand what you mean that for $r_0\le r \le r_1$ and calculated in $r_1$ , $\partial_r T(r,t)=T_1\sin\omega t$ ? – doraemonpaul Oct 12 '12 at 4:43

Of course we use separation of variables:

Let $T(r,t)=F(r)G(t)$ ,

Then $F(r)G'(t)=D\left(F''(r)G(t)+\dfrac{1}{r}F'(r)G(t)\right)$

$F(r)G'(t)=DG(t)\left(F''(r)+\dfrac{1}{r}F'(r)\right)$

$\dfrac{G'(t)}{DG(t)}=\dfrac{F''(r)+\dfrac{1}{r}F'(r)}{F(r)}=-s^2$

$\begin{cases}\dfrac{G'(t)}{G(t)}=-Ds^2\\F''(r)+\dfrac{1}{r}F'(r)+s^2F(r)=0\end{cases}$

$\begin{cases}G(t)=c_3(s)e^{-Dts^2}\\F(r)=\begin{cases}c_1(s)J_0(rs)+c_2(s)Y_0(rs)&\text{when}~s\neq0\\c_1\ln r+c_2&\text{when}~s=0\end{cases}\end{cases}$

$\therefore T(r,t)=C_1\ln r+C_2+\int_0^\infty C_3(s)e^{-Dts^2}J_0(rs)~ds+\int_0^\infty C_4(s)e^{-Dts^2}Y_0(rs)~ds$

$\partial_rT(r,t)=\dfrac{C_1}{r}-\int_0^\infty sC_3(s)e^{-Dts^2}J_1(rs)~ds-\int_0^\infty sC_4(s)e^{-Dts^2}Y_1(rs)~ds$

Substituting $T(r,0)=T_0$ and $\partial_rT(r_2,t)=T_1\sin\omega t$ respectively , you will get this system of equations $\begin{cases}C_1\ln r+C_2+\int_0^\infty C_3(s)J_0(rs)~ds+\int_0^\infty C_4(s)Y_0(rs)~ds=T_0\\\dfrac{C_1}{r_2}-\int_0^\infty sC_3(s)e^{-Dts^2}J_1(r_2s)~ds-\int_0^\infty sC_4(s)e^{-Dts^2}Y_1(r_2s)~ds=T_1\sin\omega t\end{cases}$ .

The only thing is that $C_3(s)$ and $C_4(s)$ are in fact really difficult to find.

-

Using Laplace Transform $\hat T(r,s) = \int_0^\infty e^{-st} T(r,t) dt$ on (2) we have $$\frac{d^2 \hat T}{d r^2} + \frac{1}{r} \frac{d \hat T}{d r} - \frac{s}{D} \hat T = - \frac{T_0}{D} \tag{1}$$ with the condition $$\frac{d \hat T}{d r}(r_1,s) = \frac{T_1 \omega}{s^2 + \omega^2} \tag{2}$$

As @in_wolfram_we_trust stated, you are missing a condition, and the problem is underdetermined. Surely there is something missing on the physics of the problem. That being said, taking the change of variables $r = b \rho$, eq. (1) transforms into $$\frac{d^2 \hat T}{d r^2} + \frac{1}{r} \frac{d \hat T}{d r} - \frac{b^2 s}{D} \hat T = - \frac{b^2 T_0}{D} \tag{3}$$

Letting $b = \sqrt{\frac{D}{s}}$, eq. (3) becomes $$\frac{d^2 \hat T}{d r^2} + \frac{1}{r} \frac{d \hat T}{d r} - \hat T = - \frac{T_0}{s} \tag{4}$$ which is the nonhomogeneous Modified Bessel Equation of order zero. Then eq. (1) has as solution $$\hat T(r,s) = c_1 I_0\left(\sqrt{\tfrac{s}{D}} r\right) + c_2 K_0\left(\sqrt{\tfrac{s}{D}} r\right) + \frac{T_0}{s} \tag{5}$$ A typical physical restriction would be that the temperature should remain finite which, in the case of the full cylinder would mean $c_2 = 0$, but since the cylinder does not contain $r = 0$, nothing can be said of $c_2$ (or $c_1$ for that matter). Once again, I feel the need to emphasize that the problem is missing a condition and is underdetermined.

Condition (2) means $$\frac{T_1 \omega}{s^2 + \omega} = c_1 \sqrt{\tfrac{s}{D}} I_0'\left(\sqrt{\tfrac{s}{D}} r_1\right) + c_2 \sqrt{\tfrac{s}{D}} K_0'\left(\sqrt{\tfrac{s}{D}} r_1\right),$$ and using known recurrence relations $$\frac{T_1 \omega}{s^2 + \omega} = c_1 \sqrt{\tfrac{s}{D}} I_1\left(\sqrt{\tfrac{s}{D}} r_1\right) - c_2 \sqrt{\tfrac{s}{D}} K_1\left(\sqrt{\tfrac{s}{D}} r_1\right) \tag{6}$$ From eq. (6), the solution on Laplace space is $$\hat T(r,s) = c_1 I_0\left(\sqrt{\tfrac{s}{D}} r\right) +\left[c_1 I_1\left(\sqrt{\tfrac{s}{D}} r_1\right) - \sqrt{\tfrac{D}{s}}\tfrac{T_1 \omega}{s^2 + \omega^2} \right]\tfrac{K_0\left(\sqrt{\tfrac{s}{D}} r\right)}{K_1\left(\sqrt{\tfrac{s}{D}} r_1\right)} + \frac{T_0}{s} \tag{7}$$ At this point, it's clear that no closed solution will be abailable (with this method, anyway), and given that there is a missing condition, I think this is as close to a solution as one might get. Fortunately, using Stehfest Algorithm you can have a pretty neat numerical inversion.

Doing some asymptotics one can build analytic expressions almost as good as the true solution, but given the lack of a second condition, it doesn't seem to worth the trouble (for now at least).

-