Heat equation/Solution to the 3-D Heat Equation in Cylindrical Coordinates

Definition
We are adding to the equation found in the 2-D heat equation in cylindrical coordinates, starting with the following definition: $$D:=(0,a) \times (0,b) \times (0,L)~.$$ By changing the coordinate system, we arrive at the following nonhomogeneous PDE for the heat equation: $$u_t=k \left [ \frac{1}{r} \left ( u_r + ru_{rr} \right ) + \frac{1}{r^2} u_{\theta \theta} + u_{zz} \right ] + h(r,\theta,z,t), \text{ where } (r,\theta,z) \in D, t \in (0,\infty)~.$$ We choose for the example the Robin boundary conditions and initial conditions as follows: $$ \begin{cases} \left\vert u(0,\theta,z,t) \right\vert < \infty \\ \alpha_1u(a,\theta,z,t)+\beta_1u_r(a,\theta,z,t)=0 \\ \alpha_2u(r,0,z,t)-\beta_2u_\theta(r,0,z,t)=0 \\ \alpha_3u(r,b,z,t)+\beta_3u_\theta(r,b,z,t)=0 \\ \alpha_4u(r,\theta,0,t)-\beta_4u_z(r,\theta,0,t)=0 \\ \alpha_5u(r,\theta,L,t)+\beta_5u_z(r,\theta,L,t)=0 \\ u(r,\theta,z,0)=f(r,\theta,z) \end{cases} $$

Solution
All of the boundary conditions are homogeneous, so we don't have to partition the solution into a "steady-state" portion and a "variable" portion. Otherwise, that would be the way to solve this problem.

Separate Variables
$$u(r,\theta,z,t)=R(r)\Theta(\theta)Z(z)T(t)$$ $$\Rightarrow R\Theta ZT' = k \left ( R \Theta ZT + \frac{1}{r} R' \Theta ZT + \frac{1}{r^2} R\ThetaZT + R\Theta Z''T \right )$$ $$\frac{T'}{kT} = \frac{R}{R} + \frac{1}{r} \frac{R'}{R}+\frac{1}{r^2} \frac{\Theta}{\Theta}+\frac{Z''}{Z}$$ There is a separation constant $$\gamma^2$$ that both sides of the equation are equivalent to. This yields: $$\frac{T'}{kT}=-\gamma^2 \Rightarrow {\color{Blue}T'+k\gamma^2T=0}$$ $$\frac{R}{R} + \frac{1}{r} \frac{R'}{R}+\frac{1}{r^2} \frac{\Theta}{\Theta} + \gamma^2 = -\frac{Z''}{Z} = \mu^2$$ The second equation yields the equations: $${\color{Blue}Z''+\mu^2Z=0}$$ $$\frac{R}{R}+\frac{1}{r}\frac{R'}{R}+\gamma^2-\mu^2=-\frac{1}{r^2}\frac{\Theta}{\Theta}$$ $$\Rightarrow r^2\frac{R}{R}+r\frac{R'}{R}+ \left ( \gamma^2-\mu^2 \right ) r^2=-\frac{\Theta}{\Theta}=\rho^2$$ This yields the following equations: $${\color{Blue}\Theta''+\rho^2 \Theta =0}$$ $${\color{Blue}r^2R''+rR'+ \left [ \left ( \gamma^2-\mu^2 \right ) r^2 - \rho^2 \right ] R =0}$$

Translate Boundary Conditions
Just like in the 2-D heat equation, the boundary conditions yield: $$ \begin{cases} \left\vert R(0) \right\vert < \infty \\ \alpha_1R(a)+\beta_1R'(a)=0 \\ \alpha_2\Theta(0)-\beta_2\Theta'(0)=0 \\ \alpha_3\Theta(b)+\beta_3\Theta'(b)=0 \\ \alpha_4Z(0)-\beta_4Z'(0)=0 \\ \alpha_5Z(L)+\beta_5Z'(L)=0 \\ \end{cases} $$

Solve SLPs
$$ \left. \begin{align} & Z''+\mu^2Z=0 \\ & \alpha_4Z(0)-\beta_4Z'(0)=0 \\ & \alpha_5Z(L)+\beta_5Z'(L)=0 \end{align} \right \} \begin{align} & \text{Eigenvalues } \mu_k \text{: solutions to equation } ( \alpha_4\alpha_5-\beta_4\beta_5\mu^2 ) \sin (\mu L)+(\alpha_4\beta_5+\alpha_5\beta_4)\mu \cos (\mu L) = 0 \\ & Z_k(z)=\beta_4\mu_k \cos (\mu_kz)+\alpha_4 \sin (\mu_kz), \; k=0,1,2,\cdots \end{align} $$

$$ \left. \begin{align} & \Theta''+\rho^2 \Theta =0 \\ & \alpha_2\Theta(0)-\beta_2\Theta'(0)=0 \\ & \alpha_3\Theta(b)+\beta_3\Theta'(b)=0 \end{align} \right \} \begin{align} & \text{Eigenvalues } \rho_m \text{: solutions to equation } ( \alpha_2\alpha_3-\beta_2\beta_3\rho^2 ) \sin (\rho b)+(\alpha_2\beta_3+\alpha_3\beta_2)\rho \cos (\rho b) = 0 \\ & \Theta_m(\theta)=\beta_2\rho_m \cos (\rho_m\theta)+\alpha_2 \sin (\rho_m\theta), \; m=0,1,2,\cdots \end{align} $$

$$ \left. \begin{align} & r^2R''+rR'+ \left [ \left ( \gamma^2-\mu^2 \right ) r^2 - \rho^2 \right ] R =0 \\ & \left\vert R(0) \right\vert < \infty \\ & \alpha_1R(a)+\beta_1R'(a)=0 \end{align} \right \}

\begin{align} & \text{Substitute } \lambda^2=\gamma^2-\mu_k^2 \text{ and } \nu^2=\rho_m, \; k,m=0,1,2,\cdots \\ & \text{Eigenvals } \lambda_{kmn} \text{: solns to eqn } ( \alpha_1\lambda a+\beta_1\rho_m ) J_{\rho_m} (\lambda a)-\beta_1 a \lambda J_{\rho_m+1} (\lambda a) = 0 \\ & R_{kmn}(r)=J_{\rho_m}(\lambda_{kmn}r), \; k,m,n=0,1,2,\cdots \\ & \gamma_{kmn}^2=\lambda_{kmn}^2+\mu_{k}^2 \end{align} $$

Solve Time Equation
$$T'+k\gamma_{kmn}^2T=0$$ The solution to the equation is: $$T_{kmn}(t)=C_{kmn}e^{-k \left ( \lambda_{kmn}^2+\mu_k^2 \right ) t}, \; k,m,n=0,1,2,\cdots$$

Step 2: Satisfy Initial Condition
Define: $$u(r,\theta,z,t)=\sum_{k,m,n=0}^\infty T_{kmn}(t)R_{kmn}(r)\Theta_m(\theta)Z_k(z)$$ Applying the initial condition: $$ \begin{align} u(r,\theta,z,0) & = f(r,\theta,z) \\ & = \sum_{k,m,n=0}^\infty T_{kmn}(0)R_{kmn}(r)\Theta_m(\theta)Z_k(z) \\ & = \sum_{k,m,n=0}^\infty C_{kmn}R_{kmn}(r)\Theta_m(\theta)Z_k(z) \\ \end{align} $$ This is the orthogonal expansion of $$f$$ in terms of $$R_{kmn} \otimes \Theta_m(\theta) \otimes Z_k~.$$ Hence, $$C_{kmn}=\frac{\int\limits_0^L \int\limits_0^b \int\limits_0^a f(r,\theta,z)R_{kmn}(r)\Theta_m(\theta)Z_k(z) r dr d\theta dz}{\int\limits_0^a R_{kmn}^2(r) r dr \int\limits_0^b \Theta_m^2(\theta) d\theta \int\limits_0^L Z_k^2(z) dz}, \; k,m,n=0,1,2,\cdots$$

Step 3: Solve the Non-homogeneous Equation
Let: $$u(r,\theta,z,t)=\sum_{k,m,n=0}^\infty T_{kmn}(t)R_{kmn}(r)\Theta_m(\theta)Z_k(z)$$ $$h(r,\theta,z,t)=\sum_{k,m,n=0}^\infty H_{kmn}(t)R_{kmn}(r)\Theta_m(\theta)Z_k(z), \; H_{kmn}(t)=\frac{\int\limits_0^L \int\limits_0^b \int\limits_0^a h(r,\theta,z,t)R_{kmn}(r)\Theta_m(\theta)Z_k(z) r dr d\theta dz}{\int\limits_0^a R_{kmn}^2(r) r dr \int\limits_0^b \Theta_m^2(\theta) d\theta \int\limits_0^L Z_k^2(z) dz}, \; k,m,n=0,1,2,\cdots$$ Substitute the expansions for u and h into the non-homogeneous equation: $$\sum T_{kmn}'R_{kmn}\Theta_mZ_k = k \left \{ \sum T_{kmn} \left [ \underbrace{\left ( R_{kmn} + \frac{1}{r} R_{kmn}' \right )}_{- \left [ \left ( \lambda_{kmn}^2 - \mu_k^2 \right ) - \frac{1}{r^2} \rho_m^2 \right ] R_{kmn}} \Theta_m Z_k + \frac{1}{r^2} R_{kmn} \underbrace{\Theta_m}_{-\rho_m^2 \Theta_m} Z_k + R_{kmn} \Theta_m \underbrace{Z_k''}_{-\mu_k^2Z_k} \right ] \right \} + \sum H_{kmn} R_{kmn} \Theta_m Z_k$$ $$\Leftrightarrow \sum T_{kmn}' \left [ R_{kmn} \Theta_m Z_k \right ] = \sum \left [ (-k\gamma_{kmn}^2) T_{kmn} + H_{kmn} \right ] \left [ R_{kmn} \Theta_m Z_k \right ]$$ From the linear independence of $$R_{kmn} \otimes \Theta_m \otimes Z_k~$$: $$T_{kmn}'(t)=k\gamma_{kmn}^2T_{kmn}(t) = H_{kmn}(t), \; k,m,n=0,1,2,\cdots$$ $$T_{kmn}(t)=e^{-k\gamma_{kmn}^2 t} \int\limits_0^t e^{k\gamma_{kmn}^2 s} H_{kmn}(s) ds + C_{kmn} e^{-k\gamma_{kmn}^2 t}$$ The undetermined coefficient satisfies the initial condition: $$C_{kmn}=\frac{\int\limits_0^L \int\limits_0^b \int\limits_0^a f(r,\theta,z)R_{kmn}(r)\Theta_m(\theta)Z_k(z) r dr d\theta dz}{\int\limits_0^a R_{kmn}^2(r) r dr \int\limits_0^b \Theta_m^2(\theta) d\theta \int\limits_0^L Z_k^2(z) dz}, \; k,m,n=0,1,2,\cdots$$