Reference table used in HW
\vec{\phi } | flux (class uses \vec{q}) | vector field. thermal energy per unit time per unit area. \left [ \frac{M}{T^{3}}\right ] |
\vec{\phi }\cdot \hat{n} | flux | Flux component that is outward normal to the surface \left [ \frac{M}{T^{3}}\right ] |
Q | heat source | heat energy generated per unit volume per unit time.\left [ \frac{M}{LT^{3}}\right ] |
e | thermal energy density. Scalar field. \left [ \frac{M}{LT^{2}}\right ] | |
\rho | density | mass density of material which heat flows in.\left [ \frac{M}{L^{3}}\right ] |
c | specific heat | energy to raise temp. of unit mass by one degree Kelvin. \left [ \frac{L^{2}}{T^{2}k^{o}}\right ] |
k_{0} | Thermal conductivity | Used in flux equation q=-k_{0}\nabla u, where u is temperature. \left [ \frac{ML}{T^{3}k^{0}}\right ] |
\kappa | Thermal diffusivity | Used in heat equation \frac{\partial u}{\partial t}=\kappa \nabla u+\tilde{Q}. Where \kappa =\frac{k_{0}}{\rho c}, u is temperature. |
conservation of energy | \frac{d}{dt}\int _{V}e\left ( x,t\right ) dv=\int _{A}\bar{q}\cdot \left ( -\hat{n}\right ) dA+\int _{V}Qdv. Each term has units \left [ \frac{ML^{2}}{T^{3}}\right ] | |
Fourier law | \vec{\phi }=-k_{0}\bar{\nabla }u. Relates flux to temperature gradient. | |
\nabla | Divergence operator | A vector operator. \bar{\nabla }=\left ( \frac{\partial }{\partial x},\frac{\partial }{\partial y},\frac{\partial }{\partial x}\right ) |
Fourier law is used to relate the flux to the temperature u by \phi =-k_{0}\frac{\partial u}{\partial x} for 1D or \vec{\phi }=-k_{0}\bar{\nabla }u in general.
In addition to conduction, there is convection present. This implies there is physical material mass flowing out of the control volume carrying thermal energy with it in addition to the process of conduction. Hence the flux is adjusted by this extra amount of thermal energy motion. The amount of mass that flows out of the surface per unit time per unit area is \left ( \bar{v}\rho \right ) \equiv \left [ \frac{L}{T}\frac{M}{L^{3}}\right ] =\left [ \frac{M}{T}\frac{1}{L^{2}}\right ] . Where \rho \equiv \left [ \frac{M}{L^{3}}\right ] is the mass density of the material and \bar{v}\equiv \left [ \frac{L}{T}\right ] is velocity vector of material flow at the surface.
Amount of thermal energy that \left ( \bar{v}\rho \right ) contains is given by \left ( \bar{v}\rho \right ) cu where c is the specific heat and u is the temperature. Therefore \left ( \bar{v}\rho \right ) cu is the additional flux due to convection part. Total flux becomes \begin{equation} \vec{\phi }=-k_{0}\bar{\nabla }u+\bar{v}\rho cu \tag{1} \end{equation}
\begin{align} x & =r\cos \theta \tag{1}\\ y & =r\sin \theta \tag{2} \end{align}
since r^{2}=x^{2}+y^{2} then taking derivative w.r.t. x\begin{align} 2r\frac{\partial r}{\partial x} & =2x\nonumber \\ \frac{\partial r}{\partial x} & =\frac{x}{r}\nonumber \\ & =\frac{r\cos \theta }{r}\nonumber \\ & =\cos \theta \tag{3} \end{align}
And taking derivative w.r.t. y\begin{align} 2r\frac{\partial r}{\partial y} & =2y\nonumber \\ \frac{\partial r}{\partial y} & =\frac{y}{r}\nonumber \\ & =\frac{r\sin \theta }{r}\nonumber \\ & =\sin \theta \tag{4} \end{align}
Now taking derivative w.r.t. y of (2) gives 1=\frac{\partial r}{\partial y}\sin \theta +r\frac{\partial \sin \theta }{\partial y}
Hence \fbox{$\frac{\partial \theta }{\partial y}=\frac{\cos \theta }{r}$}
Hence \fbox{$\frac{\partial \theta }{\partial x}=\frac{\sin \theta }{r}$}
By definition of unit vector\begin{align*} \hat{r} & =\frac{\bar{r}}{\left \vert r\right \vert }=\frac{\left ( \left \vert r\right \vert \cos \theta \right ) \hat{\imath }+\left ( \left \vert r\right \vert \sin \theta \right ) \hat{\jmath }}{\left \vert r\right \vert }\\ & =\cos \theta \hat{\imath }+\sin \theta \hat{\jmath } \end{align*}
To find \hat{\theta }, two relations are used.\left \Vert \hat{\theta }\right \Vert =1 by definite of unit vector. Also \hat{\theta }\cdot \hat{r}=0 since these are orthogonal vectors (basis vectors). Assuming that \hat{\theta }=c_{1}\hat{\imath }+c_{2}\hat{\jmath }, the two equations generated are\begin{align} \left \Vert \hat{\theta }\right \Vert & =1=c_{1}^{2}+c_{2}^{2}\tag{1}\\ \hat{\theta }\cdot \hat{r} & =0=\left ( \cos \theta \hat{\imath }+\sin \theta \hat{\jmath }\right ) \cdot \left ( c_{1}\hat{\imath }+c_{2}\hat{\jmath }\right ) =c_{1}\cos \theta +c_{2}\sin \theta \tag{2} \end{align}
From (2), c_{1}=\frac{-c_{2}\sin \theta }{\cos \theta }. Substituting this into (1) gives\begin{align*} 1 & =\left ( \frac{-c_{2}\sin \theta }{\cos \theta }\right ) ^{2}+c_{2}^{2}\\ & =\frac{c_{2}^{2}\sin ^{2}\theta }{\cos ^{2}\theta }+c_{2}^{2} \end{align*}
Solving for c_{2} gives\begin{align*} \cos ^{2}\theta & =c_{2}^{2}\left ( \sin ^{2}\theta +\cos ^{2}\theta \right ) \\ c_{2} & =\cos \theta \end{align*}
Since c_{2} is now known, c_{1} is found from (2)\begin{align*} 0 & =c_{1}\cos \theta +c_{2}\sin \theta \\ 0 & =c_{1}\cos \theta +\left ( \cos \theta \right ) \sin \theta \\ c_{1} & =\frac{-\left ( \cos \theta \right ) \sin \theta }{\cos \theta } \end{align*}
Hence c_{1}=-\sin \theta . Therefore \fbox{$\hat{\theta }=-\sin \theta \hat{\imath }+\cos \theta \hat{\jmath }$}
\begin{equation} \nabla =\frac{\partial }{\partial x}\hat{\imath }+\frac{\partial }{\partial x}\hat{\jmath } \tag{1} \end{equation}
Since x\equiv x\left ( r,\theta \right ) ,y\equiv y\left ( r,\theta \right ) , then\begin{align*} \frac{\partial }{\partial x} & =\frac{\partial }{\partial r}\frac{\partial r}{\partial x}+\frac{\partial }{\partial \theta }\frac{\partial \theta }{\partial x}\\ \frac{\partial }{\partial y} & =\frac{\partial }{\partial r}\frac{\partial r}{\partial y}+\frac{\partial }{\partial \theta }\frac{\partial \theta }{\partial y} \end{align*}
Equation (1) becomes \nabla =\left ( \frac{\partial }{\partial r}\frac{\partial r}{\partial x}+\frac{\partial }{\partial \theta }\frac{\partial \theta }{\partial x}\right ) \hat{\imath }+\left ( \frac{\partial }{\partial r}\frac{\partial r}{\partial y}+\frac{\partial }{\partial \theta }\frac{\partial \theta }{\partial y}\right ) \hat{\jmath }
Using result from (b), the above simplifies to \fbox{$\nabla =\hat{r}\frac{\partial }{\partial r}+\hat{\theta }\frac{1}{r}\frac{\partial }{\partial \theta }$}
\begin{align*} \bar{A} & =A_{r}\hat{r}+A_{\theta }\hat{\theta }\\ \nabla & =\hat{r}\frac{\partial }{\partial r}+\hat{\theta }\frac{1}{r}\frac{\partial }{\partial \theta } \end{align*}
Hence\begin{align} \nabla \cdot \bar{A} & =\left ( \hat{r}\frac{\partial }{\partial r}+\hat{\theta }\frac{1}{r}\frac{\partial }{\partial \theta }\right ) \cdot \left ( A_{r}\hat{r}+A_{\theta }\hat{\theta }\right ) \nonumber \\ & =\left ( \hat{r}\frac{\partial }{\partial r}\cdot A_{r}\hat{r}\right ) +\left ( \hat{r}\frac{\partial }{\partial r}\cdot A_{\theta }\hat{\theta }\right ) +\left ( \hat{\theta }\frac{1}{r}\frac{\partial }{\partial \theta }\cdot A_{r}\hat{r}\right ) +\left ( \hat{\theta }\frac{1}{r}\frac{\partial }{\partial \theta }\cdot A_{\theta }\hat{\theta }\right ) \tag{1} \end{align}
But \begin{align} \hat{r}\frac{\partial }{\partial r}\cdot A_{r}\hat{r} & =\hat{r}\frac{\partial }{\partial r}\left ( A_{r}\hat{r}\right ) \nonumber \\ & =\hat{r}\cdot \left ( \frac{\partial A_{r}}{\partial r}\hat{r}+A_{r}\frac{\partial \hat{r}}{\partial r}\right ) \nonumber \\ & =\frac{\partial A_{r}}{\partial r}\left ( \hat{r}\cdot \hat{r}\right ) +A_{r}\left ( \hat{r}\cdot \frac{\partial \hat{r}}{\partial r}\right ) \nonumber \\ & =\frac{\partial A_{r}}{\partial r}\left ( 1\right ) +A_{r}\left ( 0\right ) \nonumber \\ & =\frac{\partial A_{r}}{\partial r} \tag{2} \end{align}
And\begin{align} \hat{r}\frac{\partial }{\partial r}\cdot A_{\theta }\hat{\theta } & =\hat{r}\frac{\partial }{\partial r}\left ( A_{\theta }\hat{\theta }\right ) \nonumber \\ & =\hat{r}\cdot \left ( \frac{\partial A_{\theta }}{\partial r}\hat{\theta }+A_{\theta }\frac{\partial \hat{\theta }}{\partial r}\right ) \nonumber \\ & =\frac{\partial A_{\theta }}{\partial r}\left ( \hat{r}\cdot \hat{\theta }\right ) +A_{\theta }\left ( \hat{r}\cdot \frac{\partial \hat{\theta }}{\partial r}\right ) \nonumber \\ & =\frac{\partial A_{\theta }}{\partial r}\left ( 0\right ) +A_{\theta }\left ( 0\right ) =0 \tag{3} \end{align}
And\begin{align*} \hat{\theta }\frac{1}{r}\frac{\partial }{\partial \theta }\cdot A_{r}\hat{r} & =\hat{\theta }\frac{1}{r}\frac{\partial }{\partial \theta }\left ( A_{r}\hat{r}\right ) \\ & =\frac{1}{r}\hat{\theta }\cdot \left ( \frac{\partial A_{r}}{\partial \theta }\hat{r}+A_{r}\frac{\partial \hat{r}}{\partial \theta }\right ) \\ & =\frac{1}{r}\frac{\partial A_{r}}{\partial \theta }\left ( \hat{\theta }\cdot \hat{r}\right ) +\frac{1}{r}A_{r}\left ( \hat{\theta }\cdot \frac{\partial \hat{r}}{\partial \theta }\right ) \\ & =\frac{1}{r}\frac{\partial A_{r}}{\partial \theta }\left ( 0\right ) +\frac{1}{r}A_{r}\left ( \hat{\theta }\cdot \hat{\theta }\right ) \end{align*}
Since \frac{\partial \hat{r}}{\partial \theta }=\hat{\theta }. Therefore\begin{equation} \hat{\theta }\frac{1}{r}\frac{\partial }{\partial \theta }\cdot A_{r}\hat{r}=\frac{1}{r}A_{r} \tag{4} \end{equation}
Substituting (2,3,4,5) into (1) gives\begin{align*} \nabla \cdot \bar{A} & =\frac{\partial A_{r}}{\partial r}+0+\frac{1}{r}A_{r}+\frac{1}{r}\frac{\partial A_{\theta }}{\partial \theta }\\ & =\frac{1}{r}A_{r}+\frac{\partial A_{r}}{\partial r}+\frac{1}{r}\frac{\partial A_{\theta }}{\partial \theta } \end{align*}
Add since \frac{\partial }{\partial r}\left ( rA_{r}\right ) =A_{r}+r\frac{\partial A_{r}}{\partial r}, the above can also be written as\begin{align*} \nabla \cdot \bar{A} & =\frac{1}{r}\left ( A_{r}+r\frac{\partial A_{r}}{\partial r}\right ) +\frac{1}{r}\frac{\partial A_{\theta }}{\partial \theta }\\ & =\frac{1}{r}\frac{\partial }{\partial r}\left ( rA_{r}\right ) +\frac{1}{r}\frac{\partial A_{\theta }}{\partial \theta } \end{align*}
From part (c), it was found that \nabla =\hat{r}\frac{\partial }{\partial r}+\hat{\theta }\frac{1}{r}\frac{\partial }{\partial \theta }
Using result of part (d), which says that \nabla \cdot \bar{A}=\frac{1}{r}\frac{\partial }{\partial r}\left ( rA_{r}\right ) +\frac{1}{r}\frac{\partial A_{\theta }}{\partial \theta }, the above becomes (where now A_{r}\equiv \frac{\partial }{\partial r},A_{\theta }\equiv \frac{1}{r}\frac{\partial }{\partial \theta })\begin{align*} \nabla ^{2} & =\frac{1}{r}\frac{\partial }{\partial r}\left ( r\frac{\partial }{\partial r}\right ) +\frac{1}{r}\frac{\partial }{\partial \theta }\left ( \frac{1}{r}\frac{\partial }{\partial \theta }\right ) \\ & =\frac{1}{r}\frac{\partial }{\partial r}\left ( r\frac{\partial }{\partial r}\right ) +\frac{1}{r^{2}}\frac{\partial ^{2}}{\partial \theta ^{2}} \end{align*}
Hence \nabla ^{2}u=\frac{1}{r}\frac{\partial }{\partial r}\left ( r\frac{\partial u}{\partial r}\right ) +\frac{1}{r^{2}}\frac{\partial ^{2}u}{\partial \theta ^{2}}
Let u\equiv u\left ( r\right ) . From problem 2 part (a) it was found that\begin{align*} x & =r\cos \theta \\ y & =r\sin \theta \\ \frac{\partial r}{\partial x} & =\cos \theta \\ \frac{\partial r}{\partial y} & =\sin \theta \\ \frac{\partial \theta }{\partial y} & =\frac{\cos \theta }{r}\\ \frac{\partial \theta }{\partial x} & =\frac{-\sin \theta }{r} \end{align*}
And\begin{equation} \nabla ^{2}u=\frac{\partial ^{2}u}{\partial x^{2}}+\frac{\partial ^{2}u}{\partial y^{2}} \tag{1} \end{equation}
And\begin{align} \frac{\partial ^{2}u}{\partial y^{2}} & =\frac{\partial }{\partial y}\left ( \frac{\partial u}{\partial y}\right ) \nonumber \\ & =\frac{\partial }{\partial y}\left ( \frac{\partial u}{\partial r}\frac{\partial r}{\partial y}\right ) \nonumber \\ & =\frac{\partial }{\partial y}\left ( \frac{\partial u}{\partial r}\sin \theta \right ) \nonumber \\ & =\left ( \frac{\partial }{\partial y}\frac{\partial u}{\partial r}\right ) \sin \theta +\frac{\partial u}{\partial r}\frac{\partial \sin \theta }{\partial y}\nonumber \\ & =\left ( \frac{\partial ^{2}u}{\partial r^{2}}\frac{\partial r}{\partial y}\right ) \sin \theta +\frac{\partial u}{\partial r}\left ( \cos \theta \frac{\partial \theta }{\partial y}\right ) \nonumber \\ & =\left ( \frac{\partial ^{2}u}{\partial r^{2}}\sin \theta \right ) \sin \theta +\frac{\partial u}{\partial r}\left ( \cos \theta \left ( \frac{\cos \theta }{r}\right ) \right ) \nonumber \\ & =\frac{\partial ^{2}u}{\partial r^{2}}\sin ^{2}\theta +\frac{1}{r}\cos ^{2}\theta \frac{\partial u}{\partial r} \tag{3} \end{align}
Substituting (2),(3) into (1) gives\begin{align*} \nabla ^{2}u & =\left ( \frac{\partial ^{2}u}{\partial r^{2}}\cos ^{2}\theta +\frac{1}{r}\sin ^{2}\theta \frac{\partial u}{\partial r}\right ) +\left ( \frac{\partial ^{2}u}{\partial r^{2}}\sin ^{2}\theta +\frac{1}{r}\cos ^{2}\theta \frac{\partial u}{\partial r}\right ) \\ & =\frac{\partial ^{2}u}{\partial r^{2}}+\frac{1}{r}\left ( \sin ^{2}\theta \frac{\partial u}{\partial r}+\cos ^{2}\theta \frac{\partial u}{\partial r}\right ) \\ & =\frac{\partial ^{2}u}{\partial r^{2}}+\frac{1}{r}\frac{\partial u}{\partial r} \end{align*}
Which can be written as \fbox{$\nabla ^2u=\frac{1}{r}\frac{\partial }{\partial r}\left ( r\frac{\partial u}{\partial r}\right ) $}
Considering the thermal energy in a annulus as shown
Integrating gives total thermal energy \begin{align*} \int _{0}^{2\pi }\int _{a}^{b}\left ( c\rho u\right ) rdrd\theta & =\int _{0}^{2\pi }d\theta \int _{a}^{b}\left ( c\rho u\right ) rdr\\ & =2\pi \int _{a}^{b}\left ( c\rho u\right ) rdr \end{align*}
Using Fourier law, \begin{align*} \vec{\phi } & =-k_{0}\bar{\nabla }u\\ & =-k_{0}\left ( \hat{r}\frac{\partial u}{\partial r}+\hat{\theta }\frac{1}{r}\frac{\partial u}{\partial \theta }\right ) \end{align*}
Since symmetric in \theta , then \frac{\partial u}{\partial \theta }=0 and the above reduces to \vec{\phi }=-k_{0}\hat{r}\frac{\partial u}{\partial r}
But \hat{n}=\hat{r} since radial unit vector. The above becomes \int _{0}^{2\pi }-k_{0}\frac{\partial u}{\partial r}rd\theta =-\left ( 2\pi k_{0}\right ) r\frac{\partial u}{\partial r}
Applying that the rate of time change of total energy equal to flux through the boundaries gives\begin{align*} \frac{d}{dt}\left ( 2\pi \int _{a}^{b}\left ( c\rho u\right ) rdr\right ) & =-\left ( 2\pi k_{0}\right ) a\left . \frac{\partial u}{\partial r}\right \vert _{r=a}+\left ( 2\pi k_{0}\right ) b\left . \frac{\partial u}{\partial r}\right \vert _{r=b}\\ & =2\pi k_{0}\int _{a}^{b}\frac{\partial }{\partial r}\left ( r\frac{\partial u}{\partial r}\right ) dr \end{align*}
Moving \frac{d}{dt} inside the first integral, it become partial 2\pi \int _{a}^{b}\left ( c\rho \frac{\partial u}{\partial t}\right ) rdr=2\pi k_{0}\int _{a}^{b}\frac{\partial }{\partial r}\left ( r\frac{\partial u}{\partial r}\right ) dr
Hence \fbox{$\frac{\partial u}{\partial t}=\frac{\kappa }{r}\frac{\partial }{\partial r}\left ( r\frac{\partial u}{\partial r}\right ) $}
The earlier problem is now repeated but in this problem c\equiv c\left ( r\right ) ,k_{0}\equiv k_{0}\left ( r\right ) and \rho \equiv \rho \left ( r\right ) . These are the thermal properties in the problem.
\begin{align*} \int _{0}^{2\pi }\int _{a}^{b}\left ( c\left ( r\right ) \rho \left ( r\right ) u\right ) rdrd\theta & =\int _{0}^{2\pi }d\theta \int _{a}^{b}\left ( c\left ( r\right ) \rho \left ( r\right ) u\right ) rdr\\ & =2\pi \int _{a}^{b}\left ( c\left ( r\right ) \rho \left ( r\right ) u\right ) rdr \end{align*}
\vec{\phi }=-k_{0}\left ( r\right ) \hat{r}\frac{\partial u}{\partial r}
Applying that the rate of time change of total energy equal to flux through the boundaries gives
\begin{align*} \frac{d}{dt}\left ( 2\pi \int _{a}^{b}\left ( c\left ( r\right ) \rho \left ( r\right ) u\right ) rdr\right ) & =-\left ( 2\pi \left . k_{0}\right \vert _{r=a}\right ) a\left . \frac{\partial u}{\partial r}\right \vert _{r=a}+\left ( 2\pi \left . k_{0}\right \vert _{r=b}\right ) b\left . \frac{\partial u}{\partial r}\right \vert _{r=b}\\ & =2\pi \int _{a}^{b}\frac{\partial }{\partial r}\left ( k_{0}\left ( r\right ) r\frac{\partial u}{\partial r}\right ) dr \end{align*}
Moving \frac{d}{dt} inside the first integral, it become partial 2\pi \int _{a}^{b}\left ( c\left ( r\right ) \rho \left ( r\right ) \frac{\partial u}{\partial t}\right ) rdr=2\pi \int _{a}^{b}\frac{\partial }{\partial r}\left ( k_{0}\left ( r\right ) r\frac{\partial u}{\partial r}\right ) dr
The heat equation is \frac{\partial u}{\partial t}=\frac{\kappa }{r}\frac{\partial }{\partial r}\left ( r\frac{\partial u}{\partial r}\right ) . At steady state \frac{\partial u}{\partial t}=0. And since circular region, symmetry in \theta is assumed and therefore temperature u depends only on r only. This means u\left ( r_{0}\right ) is the same at any angle \theta for that specific r_{0}. This becomes a second order ODE\begin{align*} \frac{\kappa }{r}\frac{d}{dr}\left ( r\frac{du}{dr}\right ) & =0\\ \frac{\kappa }{r}\left ( \frac{du}{dr}+r\frac{d^{2}u}{dr^{2}}\right ) & =0\\ \frac{d^{2}u}{dr^{2}}+\frac{1}{r}\frac{du}{dr} & =0 \end{align*}
Since \frac{\kappa }{r}\neq 0. Assuming \frac{du}{dr}=v\left ( r\right ) , the above becomes\begin{align*} \frac{dv}{dr}+\frac{1}{r}v & =0\\ \frac{dv}{dr} & =-\frac{1}{r}v\\ \frac{dv}{v} & =-\frac{dr}{r} \end{align*}
Integrating \begin{align*} \ln v & =-\ln r+c_{1}\\ v & =e^{-\ln r+c_{1}}\\ & =c_{2}e^{-\ln r}\\ & =c_{2}\frac{1}{r} \end{align*}
Where c_{2}=e^{c_{1}}. Since \frac{du}{dr}=v, then\begin{align*} \frac{du}{dr} & =c_{2}\frac{1}{r}\\ du & =c_{2}\frac{1}{r}dr \end{align*}
Integrating u\left ( r\right ) =c_{2}\ln r+c_{3}
From first equation, c_{3}=T_{1}-c_{2}\ln r_{1}. Substituting in second equation gives\begin{align*} T_{2} & =c_{2}\ln r_{2}+T_{1}-c_{2}\ln r_{1}\\ & =c_{2}\left ( \ln r_{2}-\ln r_{1}\right ) +T_{1} \end{align*}
Therefore c_{2}=\frac{T_{2}-T_{1}}{\ln r_{2}-\ln r_{1}}
Hence \fbox{$u\left ( r\right ) =T_1+\left ( T_2-T_1\right ) \frac{\ln \left ( \frac{r}{r_1}\right ) }{\ln \left ( \frac{r_2}{r_1}\right ) }$}
Insulated condition implies \frac{du}{dr}=0. So the above is repeated, but this new boundary condition is now used at r_{2}. Starting from the general solution found in part (a) u\left ( r\right ) =c_{2}\ln r+c_{3}
Last problem found the solution to the heat equation in polar coordinates with symmetry in \theta to be u\left ( r\right ) =c_{2}\ln r+c_{3}
For equilibrium the total rate of heat flow at r=a should be the same as at r=b. Circumference at r=a is 2\pi a and total rate of flow at r=a is given by \beta . Hence total heat flow rate at r=a is given by \left ( 2\pi a\right ) \left . \frac{\partial u}{\partial r}\right \vert _{r=a}=2\pi a\beta
Total heat energy is, by definition \begin{equation} E=\int _{V}c\rho udv \tag{1} \end{equation}
Equation (1) becomes, where now the r limits are from a to b\begin{align*} E & =\int _{a}^{b}c\rho u\left ( 4\pi r^{2}dr\right ) \\ & =4\pi \int _{a}^{b}c\rho ur^{2}dr \end{align*}
By definition, the flux at r=b is \phi _{b}=-k_{0}\left . \frac{\partial u}{\partial r}\right \vert _{r=b}
By conservation of thermal energy\begin{align*} \frac{d}{dt}E & =-4\pi a^{2}k_{0}\left . \frac{\partial u}{\partial r}\right \vert _{r=a}+4\pi b^{2}k_{0}\left . \frac{\partial u}{\partial r}\right \vert _{r=b}\\ \frac{d}{dt}\left ( 4\pi \int _{a}^{b}c\rho ur^{2}dr\right ) & =4\pi k_{0}\int _{a}^{b}\frac{\partial }{\partial r}\left ( r^{2}\frac{\partial u}{\partial r}\right ) dr\\ \int _{a}^{b}c\rho \frac{\partial u}{\partial t}r^{2}dr & =k_{0}\int _{a}^{b}\frac{\partial }{\partial r}\left ( r^{2}\frac{\partial u}{\partial r}\right ) dr \end{align*}
Moving everything into one integral \int _{a}^{b}\left [ c\rho \frac{\partial u}{\partial t}r^{2}-k_{0}\frac{\partial }{\partial r}\left ( r^{2}\frac{\partial u}{\partial r}\right ) \right ] dr=0
Therefore \fbox{$\frac{\partial u}{\partial t}=\frac{\kappa }{r^2}\frac{\partial }{\partial r}\left ( r^2\frac{\partial u}{\partial r}\right ) $}
The heat equation is \frac{\partial u}{\partial t}=\frac{\kappa }{r^{2}}\frac{\partial }{\partial r}\left ( r^{2}\frac{\partial u}{\partial r}\right ) . For steady state \frac{\partial u}{\partial t}=0 and assuming symmetry in \theta , the heat equation becomes an ODE in r\begin{align*} \frac{\kappa }{r^{2}}\frac{d}{dr}\left ( r^{2}\frac{du}{dr}\right ) & =0\\ \frac{d}{dr}\left ( r^{2}\frac{du}{dr}\right ) & =0\\ 2r\frac{du}{dr}+r^{2}\frac{d^{2}u}{dr^{2}} & =0 \end{align*}
For r\neq 0 r\frac{d^{2}u}{dr^{2}}+2\frac{du}{dr}=0
Integrating\begin{align*} \ln v & =-2\ln r+c\\ v & =e^{-2\ln r+c}\\ & =c_{1}e^{-2\ln r}\\ & =c_{1}\frac{1}{^{r^{2}}} \end{align*}
Therefore, since \frac{du}{dr}=v\left ( r\right ) then \begin{align*} \frac{du}{dr} & =c_{1}\frac{1}{^{r^{2}}}\\ du & =c_{1}\frac{dr}{^{r^{2}}} \end{align*}
Integrating \fbox{$u\left ( r\right ) =\frac{-c_1}{r}+c_2$}
From first equation, c_{1}=c_{2}, and from second equation 80=\frac{-c_{1}}{4}+c_{1}, hence \frac{3}{4}c_{1}=80 or c_{1}=\frac{\left ( 4\right ) \left ( 80\right ) }{3}=\frac{320}{3}. Therefore, the general solution becomes u\left ( r\right ) =-\frac{320}{3}\frac{1}{r}+\frac{320}{3}