I have a problem with the solution of the heat transportation partial differential equation in cylindric coordinates. Basically I have to numerically solve the equation for a pipe through wich the water flows with laminar profile. At the beginning of the pipe the water has temperature $T$ and then flows to infinity. The outside of the pipe has temperature $T_0 < T$ and the water slowly gets cooler. I use the function

T = 
    {10^(-7)*D[t[r, z], r] + 10^(-7)*r*D[t[r, z], r, r] == 
       0.1*(1 - (r^2)/0.0001)*r*D[t[r, z], z], 
     t[r, 0] == 200, t[1, z] == 0, t[r, 100] == 10}, 
    t, {r, 0, 1}, {z, 0, 100}]

and it seems to work, but then when i plot it with

ContourPlot[T[r, z], {r, 0, 1}, {z, 0, 100}]

I get the image below. I do not understand what i am doing wrong. It seems that the program ignores the boundary condition for $z=0$ and $z=100$. Also the solution seems not be stable because loocks like some kind of noise instead of smooth lines.

enter image description here

Thank you for all the suggestions for the solution of the problem! I rewrited the code along with the implementation of the Neumann boundary condition and the change of discretisation in the following way:

T = NDSolveValue[{-Laplacian[t[r, \[Theta], z], {r, \[Theta], z}, 
    "Cylindrical"]*(10^-5) + 
 0.001*(1 - (r^2))*D[t[r, \[Theta], z], z] == 
NeumannValue[0, r == 10^-7], DirichletCondition[t[r, \[Theta], z] == 10, 
z == 10^-7 && 10^-7 <= r <= 1 && 0 <= \[Theta] <= 2*\[Pi]],DirichletCondition[t[r, \[Theta], z] == 0, 
10^-7 <= z <= 10 && r == 1 && 0 <= \[Theta] <= 2*\[Pi]]}, t, {r, 10^-7, 1}, {\[Theta], 10^-7, 2*\[Pi]}, {z, 10^-7, 10}, MaxStepSize -> 0.0001]

and used the command

Plot3D[T[r, \[Theta], z] /. \[Theta] -> \[Pi], {r, 0, 1}, {z, 0, 1},PlotRange -> All]


ContourPlot[T[r, \[Theta], z] /. \[Theta] -> \[Pi]/2, {r, 0, 1}, {z, 0, 1},MaxRecursion -> 2]

to plot the solution:

3D representation

2D representation

Now is almost the right solution, but as you can see from the 3D representation, I keep getting that small bump at the beginning and also a fluctuation at the edge that brings the temperature below 0 which is obviously wrong. I was unable to figure out if there is any other command that could be implemented in the code to give a solution without those two anomalies.

  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – Michael E2
    Commented Sep 26, 2015 at 20:12
  • 1
    $\begingroup$ You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this this meta Q&A helpful $\endgroup$
    – Michael E2
    Commented Sep 26, 2015 at 20:12
  • 2
    $\begingroup$ Shouldn't you have a NeumannCondition for $ r = 0 $? $\endgroup$
    – Pillsy
    Commented Sep 26, 2015 at 20:36

1 Answer 1


the problem is likely a result of using default discretization.

Add a dT[r,z]/dr at r=0 and then repeat the NDSolve using a length of 10 instead of 100 and you will see a plot that makes sense - might help to view in 3D (as well as contour-plot) with axis-labels.

The parameters for which you're solving should be checked for consistency with the units of the thermal diffusivity.

This boundary-condition: t[r, 100] == 10 seems questionable - it would violate conservation of energy in most situations (if not all), since it amounts to telling the equations "solve for energy balance and heat transfer over a length of tube, but be sure the condition at the exit is a uniform temperature of 10".

The outlet temperature profile is going to be dictated by the heat-transfer.

As far as MMA goes (provided the equations are correct and you're not expecting too much with default discretization MMA does) I've found the NDSolve and variants to be quite reliable.

Try to visualize something physical - a tube you can hold in your hand - so diameter something like 0.02m (20mm), length 1 meter, put in realistic properties for liquid density, thermal conductivity, specific heat (all in consistent units), then simulate that and let commonsense be your guide

  • 1
    $\begingroup$ I think it'll be better if you add the corrected code :) $\endgroup$
    – xzczd
    Commented Sep 28, 2015 at 2:10
  • $\begingroup$ @xzczd it's hard to "correct" something without knowing what OP is trying to solve. The boundary conditions are conflicting and the problem has no meaning as stated. I don't find it instructive to "just make something work" for someone. OP's question needs to be recast before anyone can know what it needed. $\endgroup$
    – Paul_A
    Commented Nov 2, 2015 at 5:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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