# Errors in Solving Heat Equation

I am have trouble with what should be a fairly simple NDSolve operation for the unsteady state heat transfer in a rod. The rod starts off at 20oC, and its surface temperature is fixed at 37oC. I am fairly new to Mathematica, so I'm not sure what is going wrong! I keep getting the error "Infinite expression 1/0. encountered."'

My code is:

 a = 0.010
Cp = 2000
k = 0.1
rho = 900
alpha = k/(rho*Cp)

Tval[r_, t_] = NDSolve[{

(* PDE - 1D Radial Heat Equation *)
(1/r)*(D[(r*D[T[r, t], r]), r]) == (1/alpha)*D[T[r, t], t],

(* Boundary Conditions *)
T[a, t] == 37,
(D[T[r, t], r] /. r -> 0) == 0,

(* Initial Condition *)
T[r, 0] == 20},

(* Define Variables *)
T[r, t], {r, 0, a}, {t, 0, 1800}]

Plot[{Tval[r,0],Tval[r,100],Tval[r,200]},{r,0,a}]


If anyone can work out why I am getting division by zero that would be much appreciated!

• I'm not sure, but I would tend to think your equation has an actual singularity at $r=0$ and that it is wrong. By rod, do you mean disk? Maybe this is useful: mathworld.wolfram.com/HeatConductionEquationDisk.html. – anderstood Apr 4 '17 at 14:50
• Replace 0 by some small number, say 10^-4, as the inner boundary in r, and NDSolve will work. Define Tval[r_, t_] = T[r, t] /. Flatten@NDSolve[..., and Plot will work too. – bbgodfrey Apr 4 '17 at 15:34
• However, T[a, t] == 37 and T[r, 0] == 20 together have the effect of assigning T[a, 0] two different values, which could be a problem. – bbgodfrey Apr 4 '17 at 15:42
• @anderstood: The ODE is correct. It does have a singular point at $r = 0$, but it is a regular singular point, which means that you can find well-behaved series solutions in the neighborhood of $r = 0$. This is a common problem that has to be dealt with when doing PDEs in curvilinear coordinates; it arises from the fact that $\vec{\nabla}r$ is not well-defined at that point. – Michael Seifert Apr 4 '17 at 17:05
• @MichaelSeifert So there is a singularity, induced by the parametrisation and not the physics, interesting. Thank you for the informative comment & link. – anderstood Apr 4 '17 at 18:09

## 1 Answer

This gives a solution when boundary and initial condition get matched with a very steep rise of the temperature at r=a and r starts from a very small value, not 0:

    Tval[r_, t_] =
T[r, t] /.
First@NDSolve[{(*PDE-1D Radial Heat Equation*)(1/
r)*(D[(r*D[T[r, t], r]), r]) == (1/alpha)*D[T[r, t],t],
T[a, t] == 37 (1 - (37 - 20)/37 E^(-8 t)), (D[T[r, t], r] /. r -> 10^-8) ==
0, T[r, 0] == 20}, T, {r, 10^-8, a}, {t, 0, 1800}]
`