Terminal Fall Velocity of a Single Spherical Particle in a Newtonian Fluid
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

<< Graphics`Graphics`

Correlation between the Newton Number (also called the Drag coefficient) and the Reynolds number. Expression depends on flow type (i.e. laminar, transition, turbulent)

(Subscript[N, e])[x_] := 
      If[x < 2, 24/x, 
        If[2 ≤ x < 500, 18.5 x^-0.6, 0.45]]

Plotting the Newton Number or the Drag coefficient versus the Reynolds number. We use log-log scale.

plt1=LogLogPlot[(Subscript[N, e])[x], {x, 10^-3, 10^5}, 
          PlotRange -> All, PlotPoints -> 200000, 
          GridLines -> Automatic, 
          PlotStyle -> {Thickness[0.01], RGBColor[1, 0, 0]}, 
          Frame -> True]

Plotting the Archimedes Number (Ar= Ne x (Re^2)) versus the Reynolds number. This dimensionless number is independent of the fall velocity and needs to be computed first in order to determine the type of flow.

plt2=LogLogPlot[(Subscript[N, e])[x] x^2, {x, 10^-3, 10^4}, 
          PlotRange -> All, PlotPoints -> 200000, 
          GridLines -> Automatic, 
          PlotStyle -> {Thickness[0.01], RGBColor[0, 1, 0]}, 
          Frame -> True]

Plotting the Newton Number divided by the Reynolds number ( = (Ne/Re)) versus the Reynolds number.

plt3=LogLogPlot[(Subscript[N, e])[x]/x, {x, 10^-3, 10^4}, 
          PlotRange -> All, PlotPoints -> 200000, 
          GridLines -> Automatic, 
          PlotStyle -> {Thickness[0.01], RGBColor[0, 0, 1]}, 
          Frame -> True]

The three plots have a common intersection point.

Show[plt3[[2]], plt2[[2]], plt1[[2]]]

Turbulent case: 
tfv is the terminal fall velocity. Expression of the tfv depends on the flow type.
Dp is particle diameter. η is the fluid's viscosity. ρ1 is the fluid's density. ρ2 is the spherical particle density. g is the gravitational acceleration. Ar is the achimedes number.

η = 1.4 10^-3; 
    Dp = 5 10^-3; ρ2 = 1260; ρ1 = 1042; g = 9.81;

(Subscript[A, r]) = 
      4/3 Dp^3 ρ1 Abs[ρ2 - ρ1] g/η^2

sol1 = FindRoot[(Subscript[A, r]) ==  (Subscript[N, e])[x] x^2, {x, 500}]

tfv[Re_] := 
      If[Re < 2, 1/(18 η) Dp^2 Abs[ρ2 - ρ1] g, 
        If[Re < 500,  
          Re η/(Dp ρ1), (Sqrt[((3 Abs[ρ2 - ρ1] Dp g)/ρ)1])]]

tfv[x /. sol1]

Laminar case: 
tfv is the terminal fall velocity. Expression of the tfv depends on the flow type. For laminar flow, tfv is given by Stokes' law.
Dp is particle diameter. η is the fluid's viscosity. ρ1 is the fluid's density. ρ2 is the spherical particle density. g is the gravitational acceleration.

η = 1.68 10^-3; 
    Dp = 1.0 10^-6; ρ2 = 916; ρ1 = 1033; g = 9.81;

(Subscript[A, r]) = 
      4/3 Dp^3 ρ1 Abs[ρ2 - ρ1] g/η^2

sol2 = FindRoot[(Subscript[A, r]) ==  (Subscript[N, e])[x] x^2, {x, 500}]

tfv[Re_] := 
      If[Re < 2, 1/(18 η) Dp^2 Abs[ρ2 - ρ1] g, 
        If[Re < 500,  
          Re η/(Dp ρ1), (Sqrt[((3 Abs[ρ2 - ρ1] Dp g)/ρ)1])]]

tfv[x /. sol2]

Transition case:
tfv is the terminal fall velocity. Expression of the tfv depends on the flow type.
Dp is particle diameter. η is the fluid's viscosity. ρ1 is the fluid's density. ρ2 is the spherical particle density. g is the gravitational acceleration.

η = 1.68 10^-2; 
    Dp = 1.0 10^-2; ρ2 = 916; ρ1 = 1033; g = 9.81;

(Subscript[A, r]) = 
      4/3 Dp^3 ρ1 Abs[ρ2 - ρ1] g/η^2

sol3 = FindRoot[(Subscript[A, r]) ==  (Subscript[N, e])[x] x^2, {x, 500}]

tfv[Re_] := 
      If[Re < 2, 1/(18 η) Dp^2 Abs[ρ2 - ρ1] g, 
        If[Re < 500,  
          Re η/(Dp ρ1), (Sqrt[((3 Abs[ρ2 - ρ1] Dp g)/ρ)1])]]

tfv[x /. sol3]