Free Convection Past a Vertical Flat Plate
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement: The governing equations are taken from Analysis of Transport Phenomena, William M. Deen, OUP, 1998.

Off[General::"spell1"]
Off[NDSolve::"ndsz"]
Off[InterpolatingFunction::"dmval"]
Off[FindRoot::"lstol"]
Off[FindRoot::"cvmit"]

Case 1: Prandtl number Pr = 1

Pr = 1;

sol[ξ_, ψ_] := 
      NDSolve[{D[F[η], {η, 3}] + 
              3 F[η] D[F[η], {η, 2}]
 - 
              2 D[F[η], {η, 1}]^2 + θ[η] == 
            0, 
D[θ[η], {η, 2}] + 
              3 Pr F[η] D[θ[η], {η, 1}] ==  
            0, 
F[0] ==  0,  F'[0] ==  
            0,  θ[0] ==  
            1,  F''[0] ==  ξ, θ'[
              0] ==  ψ}, {F[η], θ[η]}, {η, 
          0, 10}]

θEnd[ξ_?NumericQ, ψ_?NumericQ] := θ[η] /. 
          sol[ξ, ψ] /. η -> 10

DFEnd[ξ_?NumericQ, ψ_?
          NumericQ] := (D[
                F[η] /. sol[ξ, ψ], {η, 
                  1}]^2 /. η -> 
              10) + (D[
                F[η] /. sol[ξ, ψ], {η, 
                  2}]^2 /. η -> 10)

SOL1 = 
      FindRoot[{θEnd[ξ, ψ] ==  0, 
          DFEnd[ξ, ψ] ==  0}, {ξ, 
          0.5}, {ψ, -1}]

f[1] = 
      D[F[η] /. sol[ξ /. SOL1, ψ /. SOL1], {η, 
          1}]

plt1 = Plot[%, {η, 0, 10}, PlotStyle -> Hue[1/5]]

t[1] = θ[η] /. 
        sol[ξ /. SOL1, ψ /. SOL1]

plt1bis = 
      Plot[%, {η, 0, 10}, PlotStyle -> Hue[1/5]]

Case 2: Prandtl number Pr = 10

Pr = 10;

sol[ξ_, ψ_] := 
      NDSolve[{D[F[η], {η, 3}] + 
              3 F[η] D[F[η], {η, 2}]
 - 
              2 D[F[η], {η, 1}]^2 + θ[η] == 
            0, 
D[θ[η], {η, 2}] + 
              3 Pr F[η] D[θ[η], {η, 1}] ==  
            0, 
F[0] ==  0,  F'[0] ==  
            0,  θ[0] ==  
            1,  F''[0] ==  ξ, θ'[
              0] ==  ψ}, {F[η], θ[η]}, {η, 
          0, 10}]

θEnd[ξ_?NumericQ, ψ_?NumericQ] := θ[η] /. 
          sol[ξ, ψ] /. η -> 10

DFEnd[ξ_?NumericQ, ψ_?
          NumericQ] := (D[
                F[η] /. sol[ξ, ψ], {η, 
                  1}]^2 /. η -> 
              10) + (D[
                F[η] /. sol[ξ, ψ], {η, 
                  2}]^2 /. η -> 10)

SOL10 = 
      FindRoot[{θEnd[ξ, ψ] ==  0, 
          DFEnd[ξ, ψ] ==  0}, {ξ, 
          0.5}, {ψ, -1}]

f[10] = 
      D[F[η] /. sol[ξ /. SOL10, ψ /. SOL10], {η, 
          1}]

plt10 = Plot[%, {η, 0, 10}, PlotStyle -> Hue[2/5]]

t[10] = θ[η] /. 
        sol[ξ /. SOL10, ψ /. SOL10]

plt10bis = 
      Plot[%, {η, 0, 5}, PlotStyle -> Hue[2/5], 
        PlotRange -> All]

Case 3: Prandtl number Pr = 0.01

Pr = 0.01;

sol[ξ_, ψ_] := 
      NDSolve[{D[F[η], {η, 3}] + 
              3 F[η] D[F[η], {η, 2}]
 - 
              2 D[F[η], {η, 1}]^2 + θ[η] == 
            0, 
D[θ[η], {η, 2}] + 
              3 Pr F[η] D[θ[η], {η, 1}] ==  
            0, 
F[0] ==  0,  F'[0] ==  
            0,  θ[0] ==  
            1,  F''[0] ==  ξ, θ'[
              0] ==  ψ}, {F[η], θ[η]}, {η, 
          0, 100}]

θEnd[ξ_?NumericQ, ψ_?NumericQ] := θ[η] /. 
          sol[ξ, ψ] /. η -> 100

DFEnd[ξ_?NumericQ, ψ_?
          NumericQ] := (D[
                F[η] /. sol[ξ, ψ], {η, 
                  1}]^2 /. η -> 
              100) + (D[
                F[η] /. sol[ξ, ψ], {η, 
                  2}]^2 /. η -> 100)

SOL001 = 
      FindRoot[{θEnd[ξ, ψ] ==  0, 
          DFEnd[ξ, ψ] ==  0}, {ξ, 
          5}, {ψ, -10}]

f[001] = 
      D[F[η] /. sol[ξ /. SOL001, ψ /. SOL001], {η, 
          1}]

plt001 = Plot[%, {η, 0, 5}, PlotStyle -> Hue[3/5]]

t[001] = θ[η] /. 
        sol[ξ /. SOL001, ψ /. SOL001]

plt001bis = 
      Plot[%, {η, 0, 5}, PlotStyle -> Hue[3/5]]

Case 4: Prandtl number Pr = 100

Pr = 100;

sol[ξ_, ψ_] := 
      NDSolve[{D[F[η], {η, 3}] + 
              3 F[η] D[F[η], {η, 2}]
 - 
              2 D[F[η], {η, 1}]^2 + θ[η] == 
            0, 
D[θ[η], {η, 2}] + 
              3 Pr F[η] D[θ[η], {η, 1}] ==  
            0, 
F[0] ==  0,  F'[0] ==  
            0,  θ[0] ==  
            1,  F''[0] ==  ξ, θ'[
              0] ==  ψ}, {F[η], θ[η]}, {η, 
          0, 10}]

θEnd[ξ_?NumericQ, ψ_?NumericQ] := θ[η] /. 
          sol[ξ, ψ] /. η -> 10

DFEnd[ξ_?NumericQ, ψ_?
          NumericQ] := (D[
                F[η] /. sol[ξ, ψ], {η, 
                  1}]^2 /. η -> 
              10) + (D[
                F[η] /. sol[ξ, ψ], {η, 
                  2}]^2 /. η -> 10)

SOL100 = 
      FindRoot[{θEnd[ξ, ψ] ==  0, 
          DFEnd[ξ, ψ] ==  0}, {ξ, 
          0.2}, {ψ, -2}]

f[100] = 
      D[F[η] /. sol[ξ /. SOL100, ψ /. SOL100], {η, 
          1}]

plt100 = 
      Plot[%, {η, 0, 10}, PlotStyle -> Hue[4/5]]

t[100] = θ[η] /. 
        sol[ξ /. SOL100, ψ /. SOL100]

plt100bis = 
      Plot[%, {η, 0, 5}, PlotStyle -> Hue[4/5], 
        PlotRange -> All]

We reproduce figures 12-5 and 12-6 pages 500-501
of Analysis of Transport Phenomena, William M. Deen, OUP, 1998.

Velocity profiles for the free convection past an isothermal vertical plate

Show[plt1, plt10, plt001, plt100, PlotRange -> {{0, 5}, {0, 0.6}}, 
      Frame -> True]

Temperature profiles for the free convection past an isothermal vertical plate

Show[plt1bis, plt10bis, plt001bis, plt100bis, 
      PlotRange -> {{0, 5}, {0, 1}}, Frame -> True]