separation of variables to find solution to heat equation (cylinder)

Question

I'm trying to model heat flow in a cylinder using the heat equation PDE where heat flow is only radial:

$$ \frac{\partial u}{\partial t} = \frac{1}{r} \frac{\partial u}{\partial r} + \frac{\partial^2 u}{\partial r^2} $$ The I.C for $0<r<1$: $$ u(r,0) = 0 $$

The B.Cs for $t>0$: $$ \lim_{r \rightarrow 0} \frac{\partial u}{\partial r} = 0 $$ $$ u(1,t) = 1 $$ Now, we use separation of variables: Assume that the solution to the PDE is given by: $$ u(r,t) = U(r)T(t) $$ Then by rearranging, we have: $$ \frac{T'(t)}{T(t)} = \frac{1}{r} \frac{U'(r)}{U(r)} + \frac{U''(r)}{U(r)} = -\lambda = -m^2 $$ Here $m>0.$ The first ODE is given by: $$ T'(t) = -m^{2}T(t) $$ which has the following solution: $$ \sum_{n=1}^{\infty} C_{n} \exp(-m^2t) $$ The second ODE is given by: $$ U''(r) +\frac{1}{r}U'(r) +m^{2}U(r) = 0 $$ According to the literature, this is a second order ODE that would result in Bessel's functions. How do you then find the solution to this ODE and then apply the initial and boundary conditions to then find $u(r,t)$?

Answer 1

If you represent the equation in the standard form of Bessel's equation, like @Александр Пальма showed for you in the comment, then you will find a general solution $$ U(r)=c_1 J_0(m r)+c_2 Y_0(m r) $$ Then after applying two properties of Bessel's functions, 1) $$ \frac{d }{d r} J_0(m r)=-m J_1(m r),\qquad \frac{d }{d r}Y_0(m r)=-m Y_1(m r) $$ and 2) $Y_1(m r)$ has a singularity at $r=0$, you can consider the boundary condition to fix the constants.

The answer: $c_2$ is obviously zero, while $c_1=1/J_0(m)$.

Answer 2

I was just working on this same problem. @user142288's answer is on the right track, but I think it overlooks some fundamental aspects and based on the comments regarding the initial and boundary conditions, overlooking these steps can lead to some confusion. Of course, there are many ways to solve a problem, this is just what I did.

The final solution, $u(r,t)$ will be a combination of the homogeneous solution (where the boundary conditions are 0) and the particular solution

$$ u(r,t) = u_h(r,t) + u_p(r,t) $$

The homogeneous case will have boundary conditions

$$ u_h(1, t) = 0 $$

$$ \frac{\partial u_h}{\partial r} |_{r=0} = 0 $$

and initial condition $$ u_h(r, 0) = f(r) $$

where $f(r)$ is an arbitrary profile (I believe it has to be smooth).

We can solve first for the homogeneous solution using separation of variables as you pointed out

$$ u_h(r,t) = U(r) T(t) $$

which when plugged-into the original PDE (and a little rearrangement) produces

$$ \frac{T'(t)}{T(t)} = \frac{\nabla^2 U(r)}{U(r)} = -m^2 $$

where $\nabla^2$ is the Laplace operator. The solutions to these two ODEs are

$$ T(t) = a \cdot \exp(-m^2 t) $$

and in cylindrical coordinates we get

$$ U(r) = c_1 J_0(mr) + c_2 Y_0(mr) $$

where $ J_0 $ and $ Y_0$ are Bessel functions. Applying your boundary condition at $r = 0$, i.e. $ \frac{\partial u}{\partial r} |_{r=0} = 0 $ gives

$$ 0 = c_1 m J_0'(0) + c_2 m Y_0'(0) $$

$J_0'(0) = 0$, while $Y_0'(0) \neq 0$, therefore $c_2 = 0$

The second homogeneous boundary condition, i.e. $u_h(1,t) = 0$ gives the equation

$$ J_0(m) = 0 $$

Because $J_0$ is oscillatory, there are an infinite number of $m$ values that satisfy this equation. Our homogeneous solution then becomes

$$ u_h(r,t) = \sum_{i=0}^\infty C_i J_0(m_i r) \exp(-m_i^2 t) $$

The initial condition $u_h(r,0)$ states

$$ u_h(r,0) = \sum_{i=1}^\infty C_i J_0(m_i r) $$

and orthogonality allows for the constants, $C_i$ to be solved for

$$ C_i = \frac{\int_V f(r) \cdot J_0(m_i r) dV} {\int_V [J_0(m_i r)]^2 dV} $$

where $V = \pi R^2 L$, the volume of a cylinder, and $dV = 2\pi L r dr$. The full homogeneous solution is then:

$$ u_h(r,t) = \sum_{i=1}^\infty \left[ J_0(m_i r) \exp(-m_i^2 t) \ \frac{\int_V f(r) \cdot J_0(m_i r) dV} {\int_V [J_0(m_i r)]^2 dV} \right] $$

The case of the homogeneous boundary conditions can be thought of physically as a cylinder with some arbitrary temperature profile (or concentration profile) and at time $t=0$ the surface of the cylinder is held at zero temperature. The homogeneous solution solves for the evolution of this temperature profile in time and space. We can see intuitively that at infinite time, the temperature everywhere will equal the surface temperature, i.e. 0. And because of the exponential decay this is indeed what is observed, and can be stated mathematically as

$$ \lim_{t \rightarrow \infty} u_h(r, t) = 0 $$

The problem originally asked is totally analogous: instead of tracking the profile as it tends toward 0, we are solving for a profile as it tends toward 1 (because the boundary condition at the surface is $u(1,t) = 1$).

The particular solution $u_p(r,t)$ needs to satisfy the original boundary conditions:

$$ u(1, t) = 1 $$

$$ \frac{\partial u}{\partial r} |_{r=0} = 0 $$

And we can easily verify that

$$ u_p(r,t) = 1 $$ is a solution to these boundary conditions.

To make the original problem homogeneous, we had to subtract the particular solution from the full solution

$$ u_h(r,t) = u(r,t) - u_p(r,t) $$

which includes the boundary and initial conditions. Which means for the homogeneous solution, the initial condition is

$$ u_h(r,0) = u(r,0) - u_p(r,0) = 0 - 1 = -1 $$

or

$$ f(r) = -1 $$

The full solution can then be written

$$ u(r,t) = u_p(r,t) + u_h(r,t)$$

$$ u(r,t) = 1 - \sum_{i=1}^\infty \left[ J_0(m_i r) \exp(-m_i^2 t) \ \frac{\int_V J_0(m_i r) dV} {\int_V [J_0(m_i r)]^2 dV} \right] $$

Some sanity checks confirm that $$ u(r,0) = 0 $$

and

$$ \lim_{t \rightarrow \infty} u(r,t) = 1 $$