- A tall circular cylinder with radius $R$ is in temperature equilibrium with the surroundings which has the temperature $0$.
- Heat is produced uniformly from time $t=0$ with the source density $q$.
- Determine the temperature distribution over a cylindrical cross section. Disregard from any transition resitance at the surface. ( Freely translated exercise )
The heat equation together with boundary/initial conditions should look like $$\cases{\dfrac{\partial u}{\partial t} - a\Delta u = \dfrac{aq}{\lambda} & (1)\\u(R,t) = 0, \: t>0 & (2) \\u(r,0) = 0, \:0<r<R & (3)} $$ where $a$ is the diffusivity constant and $\lambda$ the heat conductivity constant. We assume cylindrical symmetry and that $u$ is bounded. Let's define the operator $\mathcal{A} = - \dfrac{1}{r}\partial_r(r \partial_r)$ with domain $\mathcal{D_A} = \{v \in C^2(0,R) \: | \: v(R) = 0, v(r) \leq C \}$ for some constant $C$. Solving the ODE $\mathcal{A}v = \lambda v$ (not same $\lambda$ as above) together with the boundary condition gives the eigenfunctions $\varphi_k = J_0(\sqrt{\lambda_k}r)$ with corresponding eigenvalues $\lambda_k = (\alpha_{0k}/R)^2$ where $\alpha_{0k}$ are the zeros to the Bessel function $J_0$, $k=1,2,...$ We now assume that $u(r,t) = \sum_{k=1}^{\infty} u_k(t) \varphi_k(r)$, which inserted into $(1)$ yields the equation $$\sum_{k=1}^{\infty} [u_k'(t)+a\dfrac{\alpha_{0k}^2}{R^2}u_k(t)]\varphi_k(r) = \dfrac{aq}{\lambda}$$ A homogenous solution is $\sum_{k=1}^{\infty} c_k e^{-at \alpha_{0k}^2/R^2}$, but since $u(r,0) = 0 \Rightarrow c_k=0, \forall k$. Now to get the particular solution we run into...
The problem
Normally, we (well at least I) would project the right side on the ON-basis $\{\phi_k \}_{k=1}^{\infty}$ but this would involve integrating the stuff like $J_0$ and $J_0^2$ which may or may not be possible. However, the answer sheet provides the answer
\begin{align} u(r,t) & = \dfrac{qR^2}{4\lambda}(1-(\dfrac{r}{R})^2-8 \sum_{n=1}^{\infty} \dfrac{J_0 (\alpha_{0n}r/R)}{\alpha_{0n}^3J_1(\alpha_{0n})} e^{-\alpha_{0n}^2ta/R^2}) \\[3mm] & = \dfrac{2qaR^2}{\lambda} \sum_{n=1}^{\infty} \dfrac{1}{\alpha_{0n}^3J_1(\alpha_{0n})}(1-e^{-\alpha_{0n}^2at/R^2})J_0(\alpha_{0n}r/R) \end{align}
which suggests there is a simple approach. Does anyone see it?
I don't see where anything simple is begin applied, just the same old tough-ass Sturm-Liouville song and dance. Let's see, the heat flow vector is $$\vec q=-k\vec\nabla T$$ and then $$\vec\nabla\cdot\vec q=-k\nabla^2T=-\rho C_p\frac{\partial T}{\partial t}+\dot q$$ So that's why we started with $$\nabla^2T-\frac{\rho C_p}k\frac{\partial T}{\partial t}=-\frac{\dot q}k$$ The first step was to take care of that inhomogeneous part, $-\frac{\dot q}k$. Since it's time- and azimuthally-independent, we try to solve $$\frac1r\frac{d}{dr}\left(r\frac{dT_p}{dr}\right)=-\frac{\dot q}k$$ Integrating, $$r\frac{dT_p}{dr}=-\frac{\dot q}{2k}r^2+C_1$$ $$T_p=-\frac{\dot q}{4k}r^2+C_1\ln r+C_2$$ Since $T_p$ should be finite as $r\rightarrow0$, $C_1=0$. Then $T_p(R)=-\frac{\dot q}{4k}R^2+C_2=0$, so $$T_p(r)=\frac{\dot q}{4k}(R^2-r^2)$$ So now we just have to solve the homogeneous equation $$\nabla^2T-\frac{\rho C_p}k\frac{\partial T}{\partial t}=0$$ And add the particular solution. Given azimuthal and axial independence we can write $T(r,t)=\rho(r)\tau(t)$ and separate variables as $$\frac1{r\rho}\frac{\partial}{\partial r}\left(r\frac{\partial\rho}{\partial r}\right)=\frac{\rho C_p}{k}\frac{\partial\tau}{\partial t}=-\lambda^2$$ The $\rho$ equation is $$r^2\frac{d^2\rho}{dr^2}+r\frac{d\rho}{dr}+\lambda^2r^2\rho=0$$ To satisfy the boundary conditions that $\rho(R)=0$ and is also finite at $r=0$ we require $\rho_k(r)=J_0\left(\frac{\lambda_kr}R\right)$ where $J_0(\lambda_k)=0$.
The $\tau$ equation reads $$\frac{d\tau_k}{dt}=-\frac{\lambda_kk}{\rho C_p}\tau_k$$ Oh dear, $k$ for both index and thermal conductivity TT. So $$\tau_k(t)=e^{-\frac{\lambda_kk}{\rho C_p}t}$$ Adding to the particular solution we have $$T(r,t)=\frac{\dot qR^2}{4k}\left\{1-\frac{r^2}{R^2}+\sum_{k=1}^{\infty}a_kJ_0\left(\frac{\lambda_kr}R\right)e^{-\frac{\lambda_kk}{\rho C_p}t}\right\}$$ Now comes the Sturm-Liouville stuff. I am using Murray R.Spiegel's Schaum's Outline on Fourier Analysis for these formulas I don't carry around in my head. The Bessel functions have normalization $$\int_0^1xJ_n^2(\lambda_kx)dx=\frac12\left[J_n^{\prime2}(\lambda_k)+\left(1-\frac{n^2}{\lambda_k^2}\right)J_n^2(\lambda_k)\right]$$ Here, $n=0$ and $J_0^{\prime}(x)=-J_1(x)$, $J_0(\lambda_k)=0$ is a boundary condition and using the orthogonality of the Bessel functions $$\int_0^RrJ_0\left(\frac{\lambda_kr}R\right)J_0\left(\frac{\lambda_jr}R\right)dr=\frac{R^2}2J_1^{2}(\lambda_k)\delta_{jk}$$ So we just have to integrate on both sides like this to pick out the Fourier coefficients $a_k$. We need the identity $$\frac d{dx}\left[x^nJ_n(x)\right]=x^nJ_{n-1}(x)$$ $$\int_0^RrJ_0\left(\frac{\lambda_kr}R\right)=\frac{R^2}{\lambda_k^2}\int_0^{\lambda_k}xJ_0\left(x\right)dx=\left.\frac{R^2}{\lambda_k^2}xJ_1(x)\right|_0^{\lambda_k}=\frac{R^2}{\lambda_k}J_1(\lambda_k)$$ $$\begin{align}\int_0^Rr\cdot r^2J_0\left(\frac{\lambda_kr}R\right)dr&= \frac{R^4}{\lambda_k^4}\int_0^{\lambda_k}x^3J_0(x)dx\\ &=\frac{R^4}{\lambda_k^4}\left\{\left.x^2\cdot xJ_1(x)\right|_0^{\lambda_k}-2\int_0^{\lambda_k}x^2J_1(x)dx\right\}\\ &=\frac{R^4}{\lambda_k^4}\left\{\lambda_k^3J_1(\lambda_k)-\left.2x^2J_2(x)\right|_0^{\lambda_k}\right\}\\ &=\frac{R^4}{\lambda_k^4}\left\{\lambda_k^3J_1(\lambda_k)-2\lambda_k^2J_2(\lambda_k)\right\}\\ &=\frac{R^4}{\lambda_k^4}\left\{\lambda_k^3J_1(\lambda_k)-2\lambda_k^2\left(\frac2{\lambda_k}J_1(\lambda_k)-J_0(\lambda_k)\right)\right\}\\ &=R^4\left(\frac1{\lambda_k}-\frac4{\lambda_k^3}\right)J_1(\lambda_k)\end{align}$$ So all this integration amounts to $$\frac{R^2}{\lambda_k}J_1(\lambda_k)-\frac1{R^2}R^4\left(\frac1{\lambda_k}-\frac4{\lambda_k^3}\right)J_1(\lambda_k)+a_k\frac{R^2}2J_1^{2}(\lambda_k)=0$$ To satisfy the boundary condition that $T(r,0)=0$. So that tells us both that $$a_k=-\frac8{\lambda_k^3J_1(\lambda_k)}$$ and that $$1-\frac{r^2}{R^2}=\sum_{k=1}^{\infty}\frac8{\lambda_k^3J_1(\lambda_k)}J_0\left(\frac{\lambda_kr}R\right)$$ So we can put together the solution $$T(r,t)=\frac{\dot qR^2}{4k}\sum_{k=1}^{\infty}\frac8{\lambda_k^3J_1(\lambda_k)}J_0\left(\frac{\lambda_kr}R\right)\left(1-e^{\frac{\lambda_kk}{\rho C_p}t}\right)$$