Solution of the natural convection boundary-layer flow near a heated flat plate

This notebook has been written in Mathematica by 

Mark J. McCready
Professor and Chair of Chemical Engineering
University of Notre Dame
Notre Dame IN 46556
USA

mjm@nd.edu
http://www.nd.edu/~mjm/


It is copyrighted to the extent allowed by whatever laws pertain to the World Wide Web and the Internet.

I would hope that as a professional courtesy, that if you use it, that this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version:  5/25/00
More recent versions of this notebook may be available at the web site:
http://www.nd.edu/~mjm/thermal_boundarylayer.nb

Reference:  F. P. Incropera and D. P. DeWitt, Fundamentals of Heat and Mass Transfer, 4th ed. Wiley, 1996.

Summary

This notebook is intended to give a first introduction to natural convection heat transfer from a solid body.  The application of interest is heat transfer from a heated sphere (a light bulb) that is enclosed in a clear plastic box.  The physical experiment is done at Notre Dame in the undergraduate chemical engineering lab (http://www.nd.edu/~chegdept/classes/cheg358.html).  The surface temperature of the light bulb is measured as a function of position around the light bulb.  These data show the development of a thermal boundary layer.

This notebook starts with the boundary layer equations for a heated, vertical flat plate under conditions of high Grashof number.  Using a similarity variable that involves the local Grashof number, the PDE's are reduced to ODE's.  The ODE's are solved for different Prandtl numbers to show the velocity and temperature profiles.

Mathematica aside

In Mathematica, it is convenient to give all expressions a "name".  I try to pick ones that are consistent with what is being done (but sometimes "temp" is used).  This assignment is done with an "=" sign.  To make an equation, a "==" is used.  This distinction is very useful in computer algebra and is employed in all of the packages with which I am familiar .

Input notation

I would enter the dynamic boundary condition from a key pad as

dc=- γ D[η[x,t],{x,2}] -
		(Subscript[ρ, 1])(-D[(Subscript[ϕ, 1])[x,y,t],t]-U1 D[(Subscript[ϕ, 1])[x,y,t],x]-g η[x,t])+
			(Subscript[ρ, 2])(-D[(Subscript[ϕ, 2])[x,y,t],t]-
					U2 D[(Subscript[ϕ, 2])[x,y,t],x]-g η[x,t])

In Standard form, which is a little shorter, you would need the typeset window to make this practical.

dc=-(γ (∂_({x,2}))η[x,t])-(Subscript[ρ, 1]) (-(∂_t)(Subscript[ϕ, 1])[x,y,t]-U1 (∂_x)(Subscript[ϕ, 1])[x,y,t]-g η[x,t])+(Subscript[ρ, 2]) (-(∂_t)(Subscript[ϕ, 2])[x,y,t]-U2 (∂_x)(Subscript[ϕ, 2])[x,y,t]-g η[x,t])

Back to start of calculations.

Finally it may be more easily read in Traditional form as shown here.  Note that this particular cell is inactive and will not be evaluated.  You can change this, if you really want to, by selecting the cell, going into the "preferences" under the Edit menu and under Cell Options, Evaluation Options, set Evaluatable to True.

dc=-(γ (((∂^2)η(x,t))/(∂(x^2))))-(Subscript[ρ, 1]) (-((∂(Subscript[ϕ, 1])[x,y,t])/(∂t))-U1 ((∂(Subscript[ϕ, 1])[x,y,t])/(∂x))-g η(x,t))+(Subscript[ρ, 2]) (-((∂(Subscript[ϕ, 2])[x,y,t])/(∂t))-U2 ((∂(Subscript[ϕ, 2])[x,y,t])/(∂x))-g η(x,t))

Back to start of calculations.

Note that you can convert any "cryptic" input expression by selecting the cell, going up to the cell menu and selecting "convert to traditional form".  This is "shift curly t" on the keypad.

Keywords

Natural Convection Heat Transfer, Thermal Boundary Layer, Boundary-Layer Flow, Similarity Variable,  High Grashof Number Heat Transfer

Problem overview

Physical situation of interest

The problem of interest is a heated sphere in a fluid of infinite extent that is otherwise not flowing.  We expect that if the sphere is heated it will be at a higher temperature than the surrounding fluid, energy will be transferred to the fluid and the resulting temperature difference within the fluid will lead to a density gradient.  If the sphere is in a uniform gravity field, you did not want pay to put your sphere on Mir, or you could not convince NASA to give you money to put it on the space shuttle, then the density gradient will lead to a body force acting on the fluid around the sphere.  The region of low density fluid will be forced upward and replaced by higher density fluid.  The question is, will the fluid flow?

Dimensional analysis for this problem would identify a Grashof number that represents the ratio of buoyancy force that would be causing flow over the viscous force that is resisting the flow.  For slightly-heated road tar, we expect a low Grashof number and thus no flow.  For a lighted light bulb in air, the Grashof number can be calculated as (10^7) or more,  this will certainly cause a strong flow.

For the case of low Gr, heat removal will be by conduction and radiation.  If you define a Nusselt number in the usual way, hD/k, the value will be two.  This is the pure condition limit.  Since the Nusselt number gives the ratio of total heat transfer to pure condition (we are not considering radiation), we expect that as the Grashof number gets large, there will be a significant fluid flow and thus significant convective enhancement.

return to conclusions

High Grashof number heat transfer

For a Grashof number that is large compared to unity, we expect that the resistance to heat transfer, alternatively the region where most of the temperature change occurs, will be confined to a thin region close to the surface.  If this region is thin compared to the macroscopic dimensions, we expect that boundary-layer approximations will give valid simplifications to the governing equations.  This will be the case.  Thus the boundary-layer equations will serve as the basic theory for the convection/conduction heat transfer from the heated sphere.

return to conclusions

Boundary-layer physics

We recall from fluid dynamics that we expect boundary-layers to occur within the fluid close to solid surfaces at high Reynolds numbers.  The current interest where the Grashof number is large is analogous.  In a momentum boundary-layer the inertia and viscous forces are about the same order of magnitude.  The governing equation thus contains the dominant inertia and viscous terms.  As the Reynolds number is increased, the boundary layer gets thinner as δ/L ~ 1/Re^1/2.  This scaling for this can be determined from the form of the equations with out solving any differential equations.  The drag on the plate can be determined to order one accuracy simply by knowing the boundary layer thickness, δ, and calculating τ = μ du/dy =~ μ U∞/δ.  

You may not be familiar with heat transfer boundary layer behavior, but a lot of similarities occur.  First now the conduction and convection terms are the same order of magnitude.  Thus convective heat transfer involves both conduction and convection.  You will notice below that the convection-conduction boundary-layer equation for energy, written in terms of T, looks quite similar to the flow direction boundary-layer equation for momentum.  The primary difference is that the energy equation does not have a term like dp/dx, unless there is a heat source within the boundary layer.  This is a source term and in the case of fluid flow is the reason the fluid is flowing (a momentum source).  

Another similarity between the present problem and fluid flow boundary layers is that we can get the basic scaling of the solution from the form of the equations.  Thus once we have the equations, can figure out how to nondimensionalize them and finally turn the PDE's into ODE's, we will know how the Nusselt number scales with Grashof number without solving the equations! 

Strictly speaking boundary-layer theory is the first order result of some clever "perturbation" expansions of the differential equations.  Leal shows this in his book, ( L. G. Leal  (1992) Laminar Flow and Convective Transport Processes,  Butterworth.)  He also shows that for a flow around a sphere, or heat transfer from a sphere, the governing equations at the lowest order (which is what we are solving) are the same for a sphere or a flat plate.  The main reason for this is because the layer is very thin compared to the radius of the sphere.  Thus "the earth is flat"! -- or at least curvature does not matter.  

If we derive the continuity, momentum and energy equations for flow past a heated flat plate we get

return to conclusions

The basic mass, momentum and energy equations for a free convection boundary - layer

Continuity

((∂u(x,y))/(∂x))+((∂v(x,y))/(∂y)) = 0

Momentum

u(x,y)((∂u(x,y))/(∂x)) + v(x,y)  ((∂u(x,y))/(∂y))=   g (T(x,y)-T0) β    +ν (((∂^2)u(x,y))/(∂(y^2)))

Energy

u(x,y)((∂T(x,y))/(∂x)) + v(x,y)((∂T(x,y))/(∂y))  =   α (((∂^2)T(x,y))/(∂(y^2)))

return to conclusions

These equations are coupled, meaning you cannot solve any of them without simultaneously solving the other two.  Further, the momentum equation is nonlinear.

Derivation of the ODE's from the original PDE's.

Do background work to set up the problem

Existence of a similarity variable.

A good reference for this section is L. G. Leal  (1992) Laminar Flow and Convective Transport Processes,  Butterworth.)  

If you encounter a fluid flow problem in a rectangular channel, you expect the characteristic length to be the channel height. For a circular pipe, it will be the diameter or radius.  For flow past a sphere, the radius is a reasonable length.  Thus we always expect to find a characteristic length, that is contained in the physics of the flow.  However, for a boundary-layer flow, where the boundary-layer is much smaller than the scale of the container, or a boundary-layer on the outside of an object where again the layer is much (say ~1/100) smaller than the device, what now is the relevant length scale?  Would nature leave us with no characteristic length?  

The answer is that there is a natural length scale.  It is just not a geometric length scale.  In the absence of a reasonable geometric length scale we look to define a "similarity" variable.  The implication is that the velocity or temperature profiles are "self-similar", and with this proper scaling can be collapsed onto the same curve.  The length scale for the normal direction within a boundary - layer is constructed from the physical properties of the flow and any other length or time variables that come to play in the problem.

For example, for the "suddenly started" plate problem, where a long and wide flat plate is started in motion, the boundary-layer thickness scales as y/(Sqrt[ν t]), where t is time and ν is the kinematic viscosity. (The boundary layer gets thicker in time as (Sqrt[ν t]).) 

To wrap this up, you always want to investigate the scaling of your problems.  If there is no geometric length scale, (i.e., problem has a infinite spatial domain), you expect the existence of a similarity variable which is a combination of the other variables and physical properties.  

Use of a similarity variable can lead to transformations that turn the PDE's into ODE's.  This can work for both linear and nonlinear equations.

The basic mass, momentum and energy equations for a free convection boundary - layer

conteq = D[u[x,y],x] +D[v[x,y],y]

umomeq = u[x,y] D[u[x,y],x] + v[x,y] D[u[x,y],y] - g β (T[x,y]- T0) - ν   D[u[x,y],{y,2}]

Teq = u[x,y] D[T[x,y],x] + v[x,y] D[T[x,y],y] - α D[T[x,y],{y,2}]

Stream function formulation to satisfy continuity

We want to turn these three coupled PDE's into 2 coupled ODE's that are nondimensionalized.  We will define a stream function from the continuity equation, define a similarity variable involving x and y.  Then transform the equations.

Note that continuity equation is satisfied automatically by a stream function formulation,
 u  = ((∂ψ[x,y])/(∂y)) ,  v  =- ((∂ψ[x,y])/(∂x))

conteq/.{u[x,y]->D[ψ[x,y],y],v[x,y]->-D[ψ[x,y],x],Derivative[a1_,a2_][u][x,y]:>  (∂_({x,a1},{y,a2}))D[ψ[x,y],y],
Derivative[a1_,a2_][v][x,y]:>  (∂_({x,a1},{y,a2}))D[-ψ[x,y],x]}

Similarity variable formulation

Now we need to do the variable transformations to make these last two equations dimensionless.  For the simplest boundary - layer problem, start up of a Couette flow, this is relatively easy.  For this problem it is a lot more complicated.

Start by defining a similarity variable. We will just choose the correct one, there are formal procedures for finding these.  (See Leal's book for how to do this.)

eta=y/x(((((Gr[x])/4)))^(1/4))

Return to conclusions

Likewise for the stream function,

psi=f[η[x,y]] 4 ν (((((Gr[x])/4)))^(1/4))

We use the Greek, in the equations, but for substitution purposes, it is convenient to have a different name, eta and psi.

Define the velocities in terms of the stream function and similarity variable

We need to use the chain rule to find all of the derivative substitutions and to find u and v. 

First find u[x,y],  u  = ((∂ψ[x,y])/(∂y)) =   ((∂ψ)/(∂η))((∂η)/(∂y))

eu1=D[psi,y]

eu1=D[psi,η[x,y]]D[η[x,y],y]

Now get the y derivative of η and substitute.

eu2 = eu1/.D[η[x,y],y]->D[eta,y]

Now we need v.
 v  =- ((∂ψ[x,y])/(∂x))  =  - ((∂ψ)/(∂η))((∂η)/(∂x))

vee1=-D[psi,x]

vee2=vee1/.D[η[x,y],x]->D[eta,x]

Now derive the 2 nonlinear ODE's that describe momentum and energy transfer

It is easiest to do the substitutions starting with the original equations. We will see if it works.

Return to conclusions

Here is the momentum equation

ueq1=umomeq/.{u[x,y]->eu2,v[x,y]->vee2,Derivative[a1_,a2_][u][x,y]:>  (∂_({x,a1},{y,a2}))eu2,
Derivative[a1_,a2_][v][x,y]:>  (∂_({x,a1},{y,a2}))vee2}

ueq2= ueq1/.{D[η[x,y],x]->D[eta,x],D[η[x,y],y]->D[eta,y],
D[η[x,y],{x,2}]->D[eta,{x,2}],D[η[x,y],{y,2}]->D[eta,{y,2}]}

ueq3=ueq2/.η[x,y]->η

ueq4 = Expand[ ueq3 x^3/ν^2]

ueq5 =  ueq4/.Gr'[x] ->3  Gr[x]/x

ueq6 = ueq5/.T[x,y]->Tstar[η] *( Ts-T0) + T0

ufinal = Simplify[ueq6/.{β-> Gr[x] ν^2/g/(Ts-T0)/x^3} ]/Gr[x]

Return to conclusions

Here is the energy equation

Teq1=Teq/.{u[x,y]->eu2,v[x,y]->vee2,Derivative[a1_,a2_][u][x,y]:>  (∂_({x,a1},{y,a2}))eu2,
Derivative[a1_,a2_][v][x,y]:>  (∂_({x,a1},{y,a2}))vee2,
Derivative[a1_,a2_][T][x,y]:>  (∂_({x,a1},{y,a2}))Tstar[η[x,y]] (Ts-T0)}

Teq2= Teq1/.{D[η[x,y],x]->D[eta,x],D[η[x,y],y]->D[eta,y],
D[η[x,y],{x,2}]->D[eta,{x,2}],D[η[x,y],{y,2}]->D[eta,{y,2}]}

Teq3=Teq2/.η[x,y]->η

Teq4 = ExpandAll[ Teq3  ]

Up to this point, we have not specified Gr[x].  To make x vanish from the equation, we need the following form for Gr[x].  Gr[x] ~ (x^3).

Teq5 =  Teq4/.Gr'[x] ->3  Gr[x]/x

Teq6=FullSimplify[ExpandAll[2 Teq5  x^2/Sqrt[Gr[x]] /ν]]

Teq7=Teq6/(T0-Ts) /.α->ν/( Pr)

Tfinal=ExpandAll[Teq7    Pr]

Return to conclusions

The final boundary-Layer ODE's are

ufinal

Tfinal

Boundary conditions

The boundary conditions are that there is no slip on the surface, and that the temp and velocity relax to 0 far away and that the temperature on the surface is fixed

The are in terms of the newest notation:

η = 0, f = f' = 0, T* = 1,
η->∞, f' = 0, T* = 0.  

Note that f is the normal velocity, f' is the tangential velocity and T* = (( (T-T0))/((Ts-T0))), where T0 is the far away temperature.

Here is a test of NDSolve without decomposing the equations.

pr=.1;

zz=NDSolve[{2 f'[η]^2 - T[η] - 3 f[η] f''[η]- f'''[η] ==0, 3 pr f[η] T'[η] + T''[η]==0,f[0]==0, f'[0]==0,T[0] ==1,f'[15]==0, T[15]==0},{f,T},η]

Plot[f'[η]/.First[zz],{η,0,15},PlotRange->All,AxesLabel->{"η","(Subscript[v, x])/(Subscript[U, 0])"}]

Plot[T[η]/.First[zz],{η,0,15},PlotRange->All,AxesLabel->{"η","T[η]"}]

pr=1;
houter =11;

zz=NDSolve[{2 f'[η]^2 - T[η] - 3 f[η] f''[η]- f'''[η] ==0, 3 pr f[η] T'[η] + T''[η]==0,f[0]==0, f'[0]==0,T[0] ==1,f'[houter]==0, T[houter]==0},{f,T},η]

zz= NDSolve[{y1'[x]==y2[x],y2'[x]==y3[x],y3'[x]==-3  y1[x] y3[x]+2  y2[x]^2 -y4[x],
y4'[x] == y5[x],y5'[x] == -3 pr (  y1[x] y5[x]),
y1[0]==0,y2[0]==0,y2[houter]==0,
y4[0]==1,y4[houter] ==0},
{y1,y2,y3,y4,y5},x];

Plot[y2[x]/.First[zz],{x,0,houter},AxesLabel->{"η","(Subscript[v, x])/(Subscript[U, 0])"}]

Plot[y4[x]/.First[zz],{x,0,houter},PlotRange->All,AxesLabel->{"η","T[η]"}]

Can we conclude anything at this point?

Return to conclusions

We have nondimensional equations for the boundary-layer flow and temperature profile that depend only on the fluid property, the Prandtl number.  Thus we expect to be able to calculate the dimensionless heat flux which is only a function of Pr.  

The dimensional heat flux, must then be given by "unscaling" the equations. Note that doing this will involve (Gr^(1/4)).  We thus conclude that this is the key scaling for convective heat transfer and this is how the data should scale (in the simplest case.)

Correlation of the experimental data

Return to conclusions

There is a recommended correlation (from Churchill) for free convection heat transfer from a sphere:

(Subscript[Nu, D]) = 2 + ((0.589 ((Subscript[Ra, D])^(1/4)))/(([1+ (((0.469/Pr))^(9/16))])^(4/9))).

(See page 504 of I &D.) 

The Nusselt and Rayleigh numbers are based on the sphere diameter.  Thus this correlation does have (Gr^(1/4))  in agreement with the scaling of the boundary-layer equations. Note that Ra≡ Gr Pr.

Numerical solution of the coupled ODE's that result after a similarity variable is introduced to the natural convection boundary - layer equations for flow near a flat plate.

Set up the problem as a system of first order ODEs.

From page 488 in I &D, the  boundary layer equations  are:

f'''+ 3 f f''- 2 ((f')^2)+  (SuperStar[T]),
(SuperStar[T])'' + 3 Pr  f (SuperStar[T])' = 0
  
These nonlinear odes do not have a known analytical solution ( for 2 particular values of Pr, a similar set of ODE's does have an analytic solution, See: L. G. Leal  (1992) Laminar flow and Convective Transport Processes, Butterworth.)  
However,they can be easily solved numerically.    Numerical solutions for ode's are often done by creating  a system of first order odes.This is done by defining
y1=f,
y2=f'
y3=f''
y4=T
y5=T',

These have to be related.  Thus we have
y1'=y2 (by definition)
y2'=y3 (by definition)
y3'=-3 y1*y3+2 (y2^2)-y4 (which is from the original ODE.)
y4'=y5
y5'=-(y1 y5)3 Pr( which is the other ODE.)

Numerical solution uses a shooting method and a Runge - Kutta Integration.  
Reference:  Numerical Recipes in Fortran, See also the notes from ChEg 258,(D. T. Leighton).  http://www.nd.edu/~dtl/cheg258/cheg258-1999/notes/l37/overheads.html

Return to conclusions

The key equations for the shooting method are:

α · δ V = -F
Vnew = Vold + δV

where α is the Jacobian Matrix obtained by varying the initial guesses for the unknown initial conditions, F is the error in the solution produced from the current initial guesses, Vold, and δV is the correction to the initial guesses for the next step. 

This turns out to be a touchy calculation.  Here are a set of initial conditions that give a solution with some trouble.  You will find that higher Pr is harder and that if youter is too large, the calculation does not work.

finit = {{.8},{-.25}};  
pr=10; (*Prandtl number*)
youter=5.5; (* outer value of η *)
deltv = .0006;(*  increment on the initial guesses used to generate the Jacobian*)

Here is the numerical solution

Initialize some variables.

finit = {{.17},{-3.4}};  (*initial condition vector for the 2 unknown derivatives, f''[0] and T'[0]*) 
pr=10; (*Prandtl number*)
youter=5; (* outer value of η *)
deltv = .0001;(* increment on the initial guesses used to generate the Jacobian*)

If you like to know how to pick out a number from a vector, look here.

finit

finit[[1,1]]

finit[[2,1]]

We need to run through the integration once to get a first values for the error.

zz=NDSolve[
{y1'[x]==y2[x],
y2'[x]==y3[x],
y3'[x]==-3  y1[x] y3[x]+2  y2[x]^2 -y4[x],
y4'[x] == y5[x],
y5'[x] == 3 pr ( - y1[x] y5[x]),y1[0]==0,y2[0]==0,y3[0]==finit[[1,1]],
y4[0]==1,y5[0] ==finit[[2,1]]},
{y1[x],y2[x],y3[x],y4[x],y5[x]},{x,0,youter}];
eps1=-(zz/.x->youter)[[1,2,2]];
eps2=-(zz/.x->youter)[[1,4,2]];
Print[ "Uerror=  ",eps1,"   f''[0]=  ",finit[[1,1]],"  Terror= ", eps2,"  Tprime[0] =   ", finit[[2,1]]];

The streamwise velocity and dimensionless temperature should be 0 as η--> ∞.  Of course we cannot not calculate out to very high values of η.  We need just high enough η.

eps1

eps2

Here is the main iteration loop.

While[Abs[eps1]+Abs[eps2]>.0001,
(*calculate the change due to changing the velocity gradient guess*)
zz=NDSolve[
{y1'[x]==y2[x],
y2'[x]==y3[x],
y3'[x]==-3  y1[x] y3[x]+2  y2[x]^2 -y4[x],
y4'[x] == y5[x],
y5'[x] == 3 pr ( - y1[x] y5[x]),y1[0]==0,y2[0]==0,y3[0]==finit[[1,1]]+deltv,
y4[0]==1,y5[0] ==finit[[2,1]]},
{y1[x],y2[x],y3[x],y4[x],y5[x]},{x,0,youter}];
eps12=-(zz/.x->youter)[[1,2,2]];
eps22=-(zz/.x->youter)[[1,4,2]];
(*calculate the change due to changing the temperature gradient guess*)
zz=NDSolve[
{y1'[x]==y2[x],
y2'[x]==y3[x],
y3'[x]==-3  y1[x] y3[x]+2  y2[x]^2 -y4[x],
y4'[x] == y5[x],
y5'[x] == 3 pr ( - y1[x] y5[x]),y1[0]==0,y2[0]==0,y3[0]==finit[[1,1]],
y4[0]==1,y5[0] ==finit[[2,1]]+deltv},
{y1[x],y2[x],y3[x],y4[x],y5[x]},{x,0,youter}];
eps13=-(zz/.x->youter)[[1,2,2]];
eps23=-(zz/.x->youter)[[1,4,2]];
(*calculate the Jacobian*)
alf={{(eps12-eps1)/deltv,(eps13-eps1)/deltv},      
          {(eps22-eps2)/deltv,(eps23-eps2)/deltv}};
(*calculate the new guesses*)
deltvnew = -Inverse[alf] . {{eps1},{eps2}};
finitnew = finit + deltvnew;
(*update everything and save the old numbers*)
finitold= finit;
eps1old = eps1;
eps2old = eps2;
finit=finitnew;
(*run the new numbers to see if the solution has converged*)
zz=NDSolve[
{y1'[x]==y2[x],
y2'[x]==y3[x],
y3'[x]==-3  y1[x] y3[x]+2  y2[x]^2 -y4[x],
y4'[x] == y5[x],
y5'[x] == 3 pr ( - y1[x] y5[x]),y1[0]==0,y2[0]==0,y3[0]==finit[[1,1]],
y4[0]==1,y5[0] ==finit[[2,1]]},
{y1[x],y2[x],y3[x],y4[x],y5[x]},{x,0,youter}];
eps1=-(zz/.x->youter)[[1,2,2]];
eps2=-(zz/.x->youter)[[1,4,2]];
Print[ "Uerror=  ",eps1,"   f''[0]=  ",finit[[1,1]],"  Terror= ", eps2,"  Tprime[0] =   ", finit[[2,1]]];]

Here are some plots of the results

Uprp1= Plot[y2[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","\!\(\*FractionBox[\(u\\\\ x\), \(2  ν\)]\) (Gr^(-(1/2)))"} ,PlotLegends->{"Pr=0.1"},PlotStyle->Purple]

Upr500= Plot[y2[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","\!\(\*FractionBox[\(u\\\\ x\), \(2  ν\)]\) (Gr^(-(1/2)))"} ,PlotLegends->{"Pr=500"},PlotStyle->Yellow]

Upr100= Plot[y2[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","\!\(\*FractionBox[\(u\\\\ x\), \(2  ν\)]\) (Gr^(-(1/2)))"} ,PlotLegends->{"Pr=100"},PlotStyle->Orange]

Upr20= Plot[y2[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","\!\(\*FractionBox[\(u\\\\ x\), \(2  ν\)]\) (Gr^(-(1/2)))"} ,PlotLegends->{"Pr=20"},PlotStyle->Green]

Upr10= Plot[y2[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","\!\(\*FractionBox[\(u\\\\ x\), \(2  ν\)]\) (Gr^(-(1/2)))"} ,PlotLegends->{"Pr=10"},PlotStyle->Red]

Upr1= Plot[y2[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","\!\(\*FractionBox[\(u\\\\ x\), \(2  ν\)]\) (Gr^(-(1/2)))"} ,PlotLegends->{"Pr=1"}]

See page 489 of I &D for a comparison.

We can plot the temperature to see how linear the profile is.  Also note where it reaches 0.  Do some plots to tell how this scales with Pr.

Tprp1=Plot[y4[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","(SuperStar[T])"}, PlotLegends->{"Pr=0.1"},PlotStyle->Purple]

Tpr500=Plot[y4[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","(SuperStar[T])"}, PlotLegends->{"Pr=500"},PlotStyle->Yellow]

Tpr100=Plot[y4[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","(SuperStar[T])"}, PlotLegends->{"Pr=100"},PlotStyle->Orange]

Tpr20=Plot[y4[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","(SuperStar[T])"}, PlotLegends->{"Pr=20"},PlotStyle->Green]

Tpr10=Plot[y4[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","(SuperStar[T])"}, PlotLegends->{"Pr=10"},PlotStyle->Red]

Tpr1=Plot[y4[x]/.zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"η","(SuperStar[T])"}, PlotLegends->{"Pr=1"}]

If we wave the plots with names (the ones that I used have values of Pr in them), then you can "Show" the results together.  (I ran this three times to make the plots with different names so that I could "Show" them.) For Pr = 0.1, 1, 10, 100.

Show[Uprp1,Upr1,Upr10, Upr20,Upr100,Upr500]

return to conclusions

Show[Tprp1,Tpr1,Tpr10,Tpr20,Tpr100,Tpr500]

Pr scaling, 1/4?  It looks good!

We can plot the temperature to see how linear the profile is.  Also note where it reaches 0.  Do some plots to tell how this scales with Pr.

anssc1000=ParametricPlot[{x  1000^(1/4), y4[x]/. zz[[1]]},{x,0,youter/(1000)^(1/4)},PlotRange->All,AxesLabel->{"η (Pr^(1/4))","(SuperStar[T])"}, PlotLegends->{"Pr = 1000"},PlotStyle->Red, AspectRatio->1]

anssc500=ParametricPlot[{x  500^(1/4), y4[x]/. zz[[1]]},{x,0,youter/(500)^(1/4)},PlotRange->All,AxesLabel->{"η","(SuperStar[T])"}, PlotLegends->{"Pr = 500"},PlotStyle->Dashed, AspectRatio->1]

anssc100=ParametricPlot[{x  100^(1/4), y4[x]/. zz[[1]]},{x,0,youter/(100)^(1/4)},PlotRange->All,AxesLabel->{"η","(SuperStar[T])"}, PlotLegends->{"Pr = 100"},PlotStyle->Purple, AspectRatio->1]

anssc10=ParametricPlot[{x  10^(1/4), y4[x]/. zz[[1]]},{x,0,youter/(10)^(1/4)},PlotRange->All,AxesLabel->{"η (Pr^(1/4))","(SuperStar[T])"}, PlotLegends->{"Pr = 10"},PlotStyle->Orange, AspectRatio->1]

Show[anssc10,anssc500,anssc100,anssc1000]

Conclusions

1.  For natural convection flows where the Grashof number is larger, a boundary layer can be expected close to a solid surface.  

2.  For situations where a natural convection boundary-layer is occurring, heat transfer will be governed by coupled energy and momentum equations. 

3.  For transport processes that occur on a semi-infinite domain, where there is no geometric length scale, it is often possible to define a (dimensionless) similarity variable that contains the natural length scale.  

4.  It is often possible to reduce PDE's to ODE's through simplifications made possible by use of the similarity variable.  

5.  From the scaling identified by the similarity variable, it is often possible to predict the macroscopic behavior without solving the differential equations.  In this case Nu ~ (Gr^(1/4)).    

6.  This prediction agrees with the recommended correlation for high Gr heat transfer.

7.  The coupled, nonlinear ODE's can be readily solved with a shooting method.

8.  Note the shape of the temperature profile and the location of the maximum velocity.