# Plotting axisymmetric solution to transient heat equation using NDSolve

I'm having some trouble plotting my results from NDSolve in an appealing way. Below is how I tried to solve the heat equation for axisymmetric, transient diffusion in the radial direction only as a function of time. There is an internal heat generation term and constant temperature at the radial boundry (ro).

ro = .3;
l = 1;
s = 25;
\[Rho] = 1050;
c = 4097;
k = .55;
q = 19305.6*Exp[.000139*t];
bc1 = T[ro, t] == s;
bc4 = (D[T[r, t], r] /. r -> 0.0001) == 0.0;
ic = T[r, 0] == s;

heat = D[T[r, t], r, r] + (1/r)*D[T[r, t], r] + q/k == (([Rho]*c)/k)*D[T[r, t], t];

sol = NDSolve[{heat, bc1, bc4, ic}, T[r, t], {r, .0001, ro}, {t, 0, 10000}]


When I plot my results, I get cartesian coordinates with radius on one axis, time on another axis, and temperature on the third axis.

Plot3D[T[r, t] /. sol, {r, .0001, .3}, {t, 0, 10000}]


This result is correct and interpretable; however, I'm looking to get a circular plot that shows me temperature at each radial position as time evolves. Any thoughts on how to achieve this?

sol = NDSolve[{heat, bc1, bc4, ic}, T, {r, .0001, ro}, {t, 0, 10000}]
f[a_, b_] := T[a, b] /. sol[[1]]
tab = Table[
ParametricPlot3D[
Evaluate[{u Cos[v], u Sin[v], f[u, j]}], {u, 0.0001, 0.3}, {v, 0,
2 Pi}, PlotRange -> {{-0.3, 0.3}, {-0.3, 0.3}, {20, 150}},
BoxRatios -> {1, 1, 1/2}, Mesh -> None], {j, 0, 10000, 100}];


tab exported as animated gif:

or color-coding temperature (rescaling 20 to 120):

tab = Table[
ParametricPlot3D[
Evaluate[{u Cos[v], u Sin[v], f[u, j]}], {u, 0.0001, 0.3}, {v, 0,
2 Pi}, PlotRange -> {{-0.3, 0.3}, {-0.3, 0.3}, {20, 150}},
BoxRatios -> {1, 1, 1/2}, Mesh -> None,
ColorFunction ->
Function[{x, y, z}, ColorData["TemperatureMap"][(z - 20)/100]],
ColorFunctionScaling -> False], {j, 0, 10000, 500}];


• This is perfect, thank you so much! One final question: is it possible to place a bar legend that tells me the temperature represented by each color? – Chase Kayrouz Mar 16 '17 at 21:35
• @ChaseKayrouz you can add a legend using Legended. Wrap the plot in Legended, e.g. Legended[plot,BarLegend[{"TemperatureMap", {20, 120}}]]. See the documentation for legend titles etc. Good luck :) – ubpdqn Mar 17 '17 at 2:48
• You can also use 'RevolutionPlot3D[f[r, 10000], {r, 0, ro}]' instead of 'ParametricPlot3D', so you don't have to write out the parameterisation explicitly. – Alexander Erlich Nov 17 '20 at 12:24

Here is sort of a quick and dirty way. We replace r with Sqrt[x^2 + y^2] and define a piece-wise function as follows:

f[rr_, tt_] = If[rr >= 0.3 || rr <= 0.0001, 25, First@T[r, t] /. sol /. {r -> rr, t -> tt}]


Using If to give the values of the function on and outside the boundaries is necessary in order that we don't evaluate the InterpolatingFunction outside of its domain. Then, we Plot3D and Manipulate over time:

Manipulate[
Plot3D[f[Sqrt[x^2 + y^2], t],
{x, -0.3, 0.3}, {y, -0.3, 0.3},
PlotRange -> {20, 140},
PerformanceGoal -> "Quality"],
{t, 0, 10000}]


• Thanks for the help, this is what I'm looking for! Unfortunately, when I try this I get "\$Aborted" in the plot window. – Chase Kayrouz Mar 16 '17 at 6:07
• @ChaseKayrouz. I see why. See the updated code. – march Mar 16 '17 at 16:04