1.4.2 Heat PDE in 2D inside a Disk, insulated boundary conditions

   1.4.2.1 Animations

Solve for \(u\left ( r,\theta ,t\right ) \), the heat PDE \[ u_{t}=k\left ( u_{rr}+\frac{1}{r}u_{r}+\frac{1}{r^{2}}u_{\theta \theta }\right ) \] where \(0<r<a\) is the radius of the disk and \(-\pi <\theta <\pi \) and \(t>0\). Initial conditions are\[ u\left ( r,\theta ,0\right ) =f\left ( r,\theta \right ) \] Boundary conditions are\begin{align*} \left . \frac{\partial u\left ( r,\theta ,t\right ) }{\partial r}\right \vert _{r=a} & =0\\ u\left ( 0,\theta ,t\right ) & <\infty \end{align*}

Boundary conditions for \(\theta \) is that it is periodic\begin{align*} u\left ( r,-\pi ,t\right ) & =u\left ( r,\pi ,t\right ) \\ u_{\theta }\left ( r,-\pi ,t\right ) & =u_{\theta }\left ( r,\pi ,t\right ) \end{align*}

Analytical solution

Same steps as above, until the solution reaches the \(R\left ( r\right ) \) Bessel ODE as in the above problem\begin{align} r^{2}R^{\prime \prime }\left ( r\right ) +rR^{\prime }\left ( r\right ) +\left ( r^{2}\lambda ^{2}-n^{2}\right ) R\left ( r\right ) & =0\qquad n=0,1,2,3,\cdots \tag{5C}\\ R^{\prime }\left ( a\right ) & =0\nonumber \\ R\left ( 0\right ) & <\infty \nonumber \end{align}

As was done in the problem above, a check is now made to see if \(\lambda =0\) is possible.

When \(n=0\)\(,\) (5C) becomes\begin{equation} r^{2}R^{\prime \prime }\left ( r\right ) +rR^{\prime }\left ( r\right ) +r^{2}\lambda ^{2}R\left ( r\right ) =0 \tag{5D} \end{equation} For \(\lambda =0\) (5D) becomes\[ r^{2}R^{\prime \prime }\left ( r\right ) +rR^{\prime }\left ( r\right ) =0 \] Whose solution is \(R\left ( r\right ) =c_{1}\ln \left ( r\right ) +c_{2}\). Since \(R\) is bounded at \(r=0\), then \(R\left ( r\right ) =c_{2}\). Now since \(R^{\prime }\left ( r\right ) =0\) at \(r=a\), then \(c_{2}\) can be any constant. Hence in this case, for \(n=0,\lambda =0\)  is an eigenvalue with constant as the eigenfunction. This was not the case in the last problem, when the boundary condition was Dirichlet \(R\left ( a\right ) =0\) and not Neumann as in this problem.

For \(\lambda \) not zero (5D) becomes\[ r^{2}R^{\prime \prime }\left ( r\right ) +rR^{\prime }\left ( r\right ) +r^{2}\lambda ^{2}R\left ( r\right ) =0 \] This is Bessel ODE of order zero. Converting it to standard form using \(t=r\lambda \) gives\[ t^{2}R^{\prime \prime }\left ( t\right ) +tR^{\prime }\left ( t\right ) +t^{2}R\left ( t\right ) =0 \] Whose solution is \(R\left ( t\right ) =c_{1}J_{0}\left ( t\right ) +c_{2}Y_{0}\left ( t\right ) \). Converting back to \(r\) gives \(R\left ( r\right ) =c_{1}J_{0}\left ( r\lambda \right ) +c_{2}Y_{0}\left ( r\lambda \right ) \). Since bounded at \(r=0\), this becomes \(R\left ( r\right ) =c_{1}J_{0}\left ( r\lambda \right ) \). B.C. \(R^{\prime }\left ( a\right ) =0\) gives \(0=c_{1}\lambda J_{0}^{\prime }\left ( \lambda a\right ) \). For non-trivial solution, \(\lambda a\) are the positive zeros of \(J_{0}^{\prime }\left ( \lambda _{m}a\right ) \) \(m=1,2,3,\cdots \) and the eigenfunction is \(R_{0}\left ( r\right ) =c_{0m}J_{0}\left ( r\lambda _{m}\right ) \) where \(\lambda _{m}a\) are zeros of \(J_{0}^{\prime }\left ( z\right ) \).  In summary, for \(n=0\), there is eigenfunction for \(\lambda =0\) which is constant. And \(\lambda \) not zero, there are \(m=1,2,3,\cdots \) eigenvalues, each has \(R_{n=0}\left ( r\right ) =c_{0m}J_{0}\left ( r\lambda _{m}\right ) \) as eigenfunctions.

Case \(n>0\). Now (5C) becomes\begin{equation} r^{2}R^{\prime \prime }\left ( r\right ) +rR^{\prime }\left ( r\right ) +\left ( r^{2}\lambda ^{2}-n^{2}\right ) R\left ( r\right ) =0\qquad n=1,2,3,\cdots \tag{5E} \end{equation} For \(\lambda =0\) (5E) becomes

\[ r^{2}R^{\prime \prime }\left ( r\right ) +rR^{\prime }\left ( r\right ) -n^{2}R\left ( r\right ) =0\qquad n=1,2,3,\cdots \] This is Euler ODE whose solution is \(R\left ( r\right ) =c_{1}r^{n}+c_{2}r^{-n}\). Since bounded at \(r=0\), then solution becomes \(R\left ( r\right ) =c_{1}r^{n}\). Now \(R^{\prime }\left ( r\right ) =c_{1}nr^{n-1}\). At \(r=a\) this becomes \(0=c_{1}na^{n-1}\) which means \(c_{1}=0\) since \(n\) is not zero. Leading to trivial solution. Hence \(\lambda =0\) is not eigenvalue in this case.

For \(\lambda \) not zero (5E) is first converted  to Bessel ODE as was done in last problem. These steps are not repeated here since they are the same. The solution is given by\[ R_{n}\left ( \lambda r\right ) =C_{n}J_{n}\left ( \lambda r\right ) +D_{n}Y_{n}\left ( \lambda r\right ) \] Since the solution is bounded at \(r=0\) and since \(Y_{n}\left ( 0\right ) \) blows up, then \(D_{n}=0\). The solution simplifies to\[ R_{n}\left ( \lambda r\right ) =C_{n}J_{n}\left ( \lambda r\right ) \] Applying the second boundary conditions \(R^{\prime }\left ( r\right ) =0\) at \(r=a\) gives\[ \left . R_{n}^{\prime }\left ( \lambda r\right ) \right \vert _{r=a}=0=\lambda C_{n}J_{n}^{\prime }\left ( \lambda a\right ) \] For non-trivial solution \(J_{n}^{\prime }\left ( \lambda a\right ) =0\) (since \(\lambda \) is not zero now). Then \(\lambda a\) are the zeros of \(J_{n}^{\prime }\left ( z\right ) \). Let the positive zeros of \(J_{n}^{\prime }\left ( z\right ) \) be \(j_{m}^{\prime }\). For \(m=1,2,3,\cdots \). Therefore \[ \lambda _{nm}=\frac{j_{m}^{\prime }}{a}\qquad n=1,2,\cdots ,m=1,2,3,\cdots \] Everything else, follows the same as above problem. The only difference is that now \(\lambda _{nm}\) obtained from the zeros of \(J_{n}^{\prime }\left ( z\right ) \), and not the zeros of \(J_{n}\left ( z\right ) \) in the last case.

In summary: When \(n=0\), then \(\lambda =0\) is possible with constant as eigenfunction for \(R\left ( r\right ) \).  In this case, the time PDE \(T_{0}^{\prime }+\lambda ^{2}kT_{0}=0\) reduces to \(T^{\prime }=0\) which implies solution \(T\left ( t\right ) =c_{1}\). When \(n=0\) but \(\lambda \) not zero, then the time ODE becomes \(T_{0m}^{\prime }+\lambda _{0m}^{2}kT_{0m}=0\) whose solution is \(T_{0m}\left ( t\right ) =e^{-\lambda _{0m}^{2}kt}\) and now there are \(m=1,2,3,\cdots \) eigenvalues, these are the zeros of \(J_{0}^{\prime }\left ( z\right ) \). When \(n>0\) then \(\lambda >0\) and the eigenvalues are the zeros of \(J_{n}^{\prime }\left ( z\right ) \) each has eigenfunction \(R_{n}\left ( r\right ) =c_{nm}J_{n}\left ( r\lambda _{nm}\right ) \).

Hence the solution is\begin{align*} u\left ( r,\theta ,t\right ) & =\overset{n=0,\lambda =0}{\overbrace{T_{0}\left ( t\right ) R_{0}\left ( r\right ) \Theta _{0}\left ( \theta \right ) }}+\overset{n=0,\lambda \neq 0}{\overbrace{T_{0m}\left ( t\right ) R_{0m}\left ( r\right ) \Theta _{0}\left ( \theta \right ) }}\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}e^{-k\lambda _{nm}^{2}t}J_{n}\left ( \lambda _{nm}r\right ) \cos n\theta +B_{nm}e^{-k\lambda _{nm}^{2}t}J_{n}\left ( \lambda _{nm}r\right ) \sin n\theta \\ & =C_{0}+\sum _{m=1}^{\infty }\overset{n=0,\lambda \neq 0}{\overbrace{A_{0m}e^{-k\lambda _{0m}^{2}t}J_{0}\left ( \lambda _{0m}r\right ) }}+\sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}e^{-k\lambda _{nm}^{2}t}J_{n}\left ( \lambda _{nm}r\right ) \cos n\theta +B_{nm}e^{-k\lambda _{nm}^{2}t}J_{n}\left ( \lambda _{nm}r\right ) \sin n\theta \end{align*}

Where the \(C_{0}\) constant has all constants for case \(n=0,\lambda =0\) combined, since each solution is constant in this case. To find \(C_{0}\), then from initial conditions and for \(n=0,\lambda =0\) only, the above becomes\[ f\left ( r,\theta \right ) =C_{0}\] Applying orthogonality noting that eigenfunction is constant in this case gives\begin{align*} \int _{-\pi }^{\pi }\int _{0}^{a}rf\left ( r,\theta \right ) drd\theta & =C_{0}\int _{-\pi }^{\pi }\int _{0}^{a}rdrd\theta \\ & =C_{0}\pi a^{2} \end{align*}

Where \(\pi a^{2}\) is the disk area.\[ C_{0}=\frac{1}{\pi a^{2}}\int _{-\pi }^{\pi }\int _{0}^{a}rf\left ( r,\theta \right ) drd\theta \] Hence the final solution is\begin{align*} u\left ( r,\theta ,t\right ) & =\overset{n=0,\lambda =0}{\overbrace{\frac{1}{\pi a^{2}}\int _{-\pi }^{\pi }\int _{0}^{a}rf\left ( r,\theta \right ) drd\theta \theta }}+\\ & \sum _{m=1}^{\infty }\overset{n=0,\lambda >0}{\overbrace{A_{0m}e^{-k\lambda _{0m}^{2}t}J_{0}\left ( \lambda _{0m}r\right ) }}+\\ & \sum _{n=1}^{\infty }\sum _{m=1}^{\infty }A_{nm}e^{-k\lambda _{nm}^{2}t}J_{n}\left ( \lambda _{nm}r\right ) \cos n\theta +B_{nm}e^{-k\lambda _{nm}^{2}t}J_{n}\left ( \lambda _{nm}r\right ) \sin n\theta \end{align*}

To find \(A_{0m}\), for \(n=0\) and \(\lambda >0\) the above becomes\[ u\left ( r,\theta ,t\right ) =\sum _{m=1}^{\infty }A_{0m}e^{-k\lambda _{0m}^{2}t}J_{0}\left ( \lambda _{0m}r\right ) \] Applying orthogonality twice, using initial conditions at \(t=0,\) then the above becomes (notice, there is implicit \(n=0\) here, which means \(\cos \left ( n\theta \right ) \) when \(n=0\)), and this is where the \(2\pi \) comes in below)\begin{align*} \int _{-\pi }^{\pi }\int _{0}^{a}rf\left ( r,\theta \right ) J_{0}\left ( \lambda _{0m}r\right ) drd\theta & =2\pi A_{0m}\int _{0}^{a}rJ_{0}^{2}\left ( \lambda _{0m}r\right ) dr\\ A_{0m} & =\frac{\int _{-\pi }^{\pi }\int _{0}^{a}rf\left ( r,\theta \right ) J_{0}\left ( \lambda _{0m}r\right ) drd\theta }{2\pi \int _{0}^{a}rJ_{0}^{2}\left ( \lambda _{0m}r\right ) dr} \end{align*}

The other constants \(A_{nm},B_{nm}\) are found as in the last problem with \(A_{nm}\) given by (6AA,6AB) and \(B_{nm}\) by (6B). Here they are again

\begin{align} \int _{0}^{a}\left ( \int _{-\pi }^{\pi }f\left ( r,\theta \right ) \cos n\theta d\theta \right ) rJ_{n}\left ( \lambda _{nm}r\right ) dr & =\pi A_{nm}\int _{0}^{a}rJ_{n}^{2}\left ( \lambda _{nm}r\right ) dr\nonumber \\ A_{nm} & =\frac{\int _{0}^{a}\left ( \int _{-\pi }^{\pi }f\left ( r,\theta \right ) \cos \left ( n\theta \right ) d\theta \right ) rJ_{n}\left ( \lambda _{nm}r\right ) dr}{\pi \int _{0}^{a}rJ_{n}^{2}\left ( \lambda _{nm}r\right ) dr}\qquad n=1,2,\cdots ,m=1,2,3,\cdots \tag{6AB} \end{align}

Similarly, \begin{equation} B_{nm}=\frac{\int _{0}^{a}\left ( \int _{-\pi }^{\pi }f\left ( r,\theta \right ) \sin \left ( n\theta \right ) d\theta \right ) rJ_{n}\left ( \lambda _{nm}r\right ) dr}{\pi \int _{0}^{a}rJ_{n}^{2}\left ( \lambda _{nm}r\right ) dr}\qquad n=1,2,\cdots ,m=1,2,3,\cdots \tag{6B} \end{equation} and \(\lambda _{nm}=\frac{j_{m}^{\prime }}{a}\) where \(j_{m}^{\prime }\) is the \(m^{th}\) zero of \(J_{n}^{\prime }\left ( z\right ) \).

1.4.2.1 Animations

Example 1

Since the disk is insulated, then total heat energy remains the same. The final heat distribution will average out over the whole disk. As time passes, heat will level out to be the same over the whole disk reaching temperature of \(0.5\) degrees everywhere. This animation is run for 1.5 seconds, using \(k=1,a=1,f\left ( r,\theta \right ) =a-r^{2}\). It used 13 terms in the series.

Example 2

This animation uses \(a=1,f\left ( r,\theta \right ) =\left ( r-r^{2}\right ) \sin \theta ,k=1\). Same number of terms above the above animation. Animation time is 0.15 seconds.

Example 3

This animation uses \(a=1,f\left ( r,\theta \right ) =\left ( r-r^{2}\right ) \sin \theta \cos \theta ,k=1\). Same number of terms above the above animation. Animation time is 0.15 seconds.