# Heat convection differential equations from 1952 - Mathematica “fails to converge”

I am trying to solve a fundamental problem in analytical convective heat transfer: laminar free convection flow and heat transfer from a flat plate parallel to the direction of the generating body force.

Brief History of the problem

Effectively: a flat plate is vertical and parallel to the direction of gravity vector. The plate is hot and the ambient is not. Heat transfer occurs from the plate to the ambient through natural convection due to density stratification.

Simon Ostrach, a distinguished scientist in the field of microgravity science solved this problem through a coupled set of equations. In Ostrach's work, these equations were solved by an IBM Card Programmed Electronic Calculator

$$F''' + 3 FF'' - 2 (F')^2 + H = 0$$ $$H'' + 2 \text{Pr} F H' = 0$$

The Boundary conditions are: $$F'(0) = F(0) = 0$$ $$H(0) = 1$$ $$F'(\infty) = H(\infty) = 0$$

Here, $F$ provides the hydrodynamic solution while $H$ provides the thermal solution with Pr being the Prandtl number which is a property of the fluid that the plate is "immersed" in.

My Mathematica code ... it runs selectively

Clear[max, Pr, T, f, η, p];
max = 50;
Pr = 0.72;
pohl = NDSolve[{f'''[η] + 3 f[η] f''[η] -
2 (f'[η])^2 + T[η] == 0,
T''[η] + 2 Pr f[η] T'[η] == 0, f[0] == f'[0] == 0,
f'[max] == 0, T[0] == 1, T[max] == 0}, {f, T}, {η, max}]

p4 = Plot[{Evaluate[f'[η] /. pohl]}, {η, 0, max},
PlotRange -> All,
PlotLabel ->
Style[Framed["Hydrodynamic development is depicted on this plot"],
10, Blue, Background -> Lighter[Yellow]], ImageSize -> Large,
BaseStyle -> {FontWeight -> "Bold", FontSize -> 18},
AxesLabel -> {"η", "f'[η]"}, PlotLegends -> "Expressions"]


For a Prandtl number of 0.72 (Air) I get a velocity profile ($F'$) as suggested by the Ostrach in his pivotal report. However, for, many Prandtl numbers, the following warning message is sometimes flashed and I get incorrect velocity profiles (negative velocities) per the publication. For instance try Pr=6.

FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations. >>

NDSolve::berr: The scaled boundary value residual error of 2.9035865095898766*^7 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.

I have experimented with the LSODA method because this system of diff eqs is stiff and LSODA has proven to be a 'magic wand' in the past. What gives? How do I select a method for this problem? I wonder if this is a problem with the method of choice (or default method with no options) or my definition of the "free stream limit" $\infty$.

Pr=0.01

Pr=0.72

Pr=0.6 (what went wrong? Warning message was displayed too...)

• I just found this. I'll experiment with this and add more detail to this problem. – dearN Jan 15 '16 at 20:18
• Mathematica v9 returns your third (erroneous) plot for Pr = 0.72 also. – Dr. belisarius Jan 15 '16 at 21:52
• Approximating max with 50 seems to be too large for Pr = 0.6 (and Pr = 0.72 in v9), max = 20 solves the problem. I guess max = 50 can also be used if one manually sets the initial guess of "Shooting" method carefully. BTW, I observed similar problem when solving Young-Laplace equation and found using the asymptotic solution at far field as the boundary a good solution. Does this set of equation owns a analytic asymptotic solution at far field? – xzczd Jan 16 '16 at 4:02
• @xzczd Yes, this (both F and H) will be asymptotic in the far field. I found a mathematica file that uses the Shooting method by expressing one of these coupled variables as an initial condition. Alas, I do not have the time to express this here because of a lack of time today! I'll try to get to it in the evening. – dearN Jan 16 '16 at 13:04
• @xzczd Could you share a link (if on SE, i.e) to the issue you were having with the Young-Laplace equation? I would be interested to look at it. I am into analytical fluid dynamics and heat transfer and mathematica has revolutionized it for me! – dearN Jan 16 '16 at 13:11

The problem is with the default starting initial conditions used by the shooting method in NDSolve. The shooting method is where FindRoot is being used internally, so the OP's error message is a strong hint that this is the problem. Getting convergence in a nonlinear system can depend greatly on the starting conditions.

Having luckily solved the system for Pr = 0.72, we can use its initial conditions as starting values for Pr = 0.6. We hope that it will be suitably close. (If not, we could have tried solving for, say, Pr = 0.66 and edged our way bit by bit to 0.6, hoping that the dependence on Pr is continuous.)

Pr = 0.72;
pohl72 =
NDSolve[{f'''[η] + 3 f[η] f''[η] - 2 (f'[η])^2 +
T[η] == 0, T''[η] + 2 Pr f[η] T'[η] == 0,
f[0] == f'[0] == 0, f'[max] == 0, T[0] == 1, T[max] == 0},
{f, T}, {η, max}];

Pr = 0.6;
pohl = NDSolve[{f'''[η] + 3 f[η] f''[η] - 2 (f'[η])^2 + T[η] == 0,
T''[η] + 2 Pr f[η] T'[η] == 0, f[0] == f'[0] == 0,
f'[max] == 0, T[0] == 1, T[max] == 0},
{f, T}, {η, max},
Method -> {"Shooting",
"StartingInitialConditions" ->
Thread[{f[0], f'[0], f''[0], T[0], T'[0]} ==
({f[0], f'[0], f''[0], T[0], T'[0]} /. First@pohl72)]}]


Plot:

Plot[{Evaluate[f'[η] /. pohl]}, {η, 0, max},
PlotRange -> All,
PlotLabel ->
Style[Framed["Hydrodynamic development is depicted on this plot"],
10, Blue, Background -> Lighter[Yellow]], ImageSize -> Large,
BaseStyle -> {FontWeight -> "Bold", FontSize -> 18},
AxesLabel -> {"η", "f'[η]"}, PlotLegends -> "Expressions"]


The problem with your system of equation is that the "zero" far field solution is unstable, hence extremely sensitive to the initial conditions. Posing the problem as an initial value problem, with the "known" conditions from the successful solution:

max = 50;
Pr = .72;
a = 0.7172594734816521
b = -0.4344414944896132
pohl = NDSolve[{f'''[\[Eta]] + 3 f[\[Eta]] f''[\[Eta]] -
2 (f'[\[Eta]])^2 + T[\[Eta]] == 0,
T''[\[Eta]] + 2 Pr f[\[Eta]] T'[\[Eta]] == 0, f[0] == f'[0] == 0,
f''[0] == a, T[0] == 1, T'[0] == b}, {f, T}, {\[Eta], max}]


now change the T' initial value by 1/2%:

b=.995b


and you see you get the essential solution but it eventually blows up:

so the shooting method needs an exceptionally good initial guess to work.

• What about taking the linear parts of these equations and using it's solutions as an initial guess? – zhk Jun 21 '16 at 15:20
• As a practical matter Michael E2's approach is what you need to do, start with a known solution and incrementally change the parameter. You might improve that by adding a critera to stop if f,f' and T,T` are sufficiently small, since you know constant zero solves the system. – george2079 Jun 21 '16 at 15:46