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}]


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

  • $\begingroup$ 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. $\endgroup$ – anderstood Apr 4 '17 at 14:50
  • 1
    $\begingroup$ 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. $\endgroup$ – bbgodfrey Apr 4 '17 at 15:34
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey Apr 4 '17 at 15:42
  • $\begingroup$ @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. $\endgroup$ – Michael Seifert Apr 4 '17 at 17:05
  • $\begingroup$ @MichaelSeifert So there is a singularity, induced by the parametrisation and not the physics, interesting. Thank you for the informative comment & link. $\endgroup$ – anderstood Apr 4 '17 at 18:09

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}]
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.