Viscous Flow

Real fluid behavior - friction, boundary layers, and turbulence

Introduction

Real fluids exhibit viscosity - internal resistance to deformation due to molecular momentum transfer. Viscous effects cause energy dissipation, velocity gradients near surfaces (boundary layers), and fundamentally different flow behavior than ideal fluids.

The Navier-Stokes equations, which include viscous terms, are the complete governing equations for Newtonian fluid flow - one of the most important and challenging equations in physics.

Viscosity

Newton's Law of Viscosity

For a Newtonian fluid, shear stress is proportional to velocity gradient:

$$\tau = \mu \frac{du}{dy}$$

μ = dynamic viscosity [Pa·s = kg/(m·s)]. This linear relationship defines Newtonian fluids.

Dynamic Viscosity (μ)

Absolute measure of fluid's resistance to shear. Water: ~0.001 Pa·s, Air: ~1.8×10⁻⁵ Pa·s

Kinematic Viscosity (ν)

ν = μ/ρ [m²/s]. Ratio of viscous to inertial effects. Appears in Reynolds number.

Non-Newtonian Fluids

Shear-thinning

Viscosity decreases with shear rate (paint, blood, ketchup)

Shear-thickening

Viscosity increases with shear rate (cornstarch in water)

Bingham Plastic

Requires yield stress to flow (toothpaste, mayonnaise)

Navier-Stokes Equations

The complete equations of motion for incompressible Newtonian fluids:

$$\rho\left(\frac{\partial \vec{v}}{\partial t} + (\vec{v} \cdot \nabla)\vec{v}\right) = -\nabla p + \mu \nabla^2 \vec{v} + \rho \vec{g}$$

$$\nabla \cdot \vec{v} = 0 \quad \text{(incompressibility)}$$

Component Form (Cartesian)

$$\rho\left(\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z}\right) = -\frac{\partial p}{\partial x} + \mu\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}\right) + \rho g_x$$

Similar equations for v and w components. These coupled, nonlinear PDEs are notoriously difficult to solve - proving existence and smoothness of solutions is a Millennium Prize Problem!

Why So Difficult?

• Nonlinearity: The convective term (v·∇)v couples velocity components
• Coupling: Pressure and velocity are coupled through continuity
• Turbulence: At high Re, solutions become chaotic with wide range of scales

Reynolds Number and Flow Regimes

Reynolds Number

The ratio of inertial to viscous forces - the key dimensionless parameter:

$$Re = \frac{\rho v L}{\mu} = \frac{v L}{\nu} = \frac{\text{Inertial forces}}{\text{Viscous forces}}$$

v = characteristic velocity, L = characteristic length, ν = kinematic viscosity

Flow Regimes

Laminar Flow

Smooth, orderly layers sliding past each other

• Pipe flow: Re < 2300
• Flat plate: Re_x < 5×10⁵
• Viscous forces dominate
• Predictable, steady (often)

Turbulent Flow

Chaotic, irregular motion with eddies

• Pipe flow: Re > 4000
• Flat plate: Re_x > 5×10⁵
• Inertial forces dominate
• Enhanced mixing, higher friction

Transition Region

2300 < Re < 4000 for pipes - flow can be either laminar or turbulent depending on disturbances. Hysteresis effects possible.

Pipe Flow

Laminar Pipe Flow (Hagen-Poiseuille)

For fully developed laminar flow in a circular pipe:

$$u(r) = \frac{R^2}{4\mu}\left(-\frac{dp}{dx}\right)\left(1 - \frac{r^2}{R^2}\right) = u_{max}\left(1 - \frac{r^2}{R^2}\right)$$

Parabolic velocity profile. Maximum velocity at centerline is twice the average.

$$Q = \frac{\pi R^4}{8\mu}\left(-\frac{dp}{dx}\right) = \frac{\pi R^4 \Delta p}{8\mu L}$$

Flow rate proportional to R⁴ - halving diameter reduces flow 16×!

Friction Factor

Dimensionless pressure drop (Darcy-Weisbach equation):

$$\Delta p = f \frac{L}{D} \frac{\rho v^2}{2} \quad \text{where} \quad h_f = f\frac{L}{D}\frac{v^2}{2g}$$

Laminar

$$f = \frac{64}{Re}$$

Turbulent (Smooth)

$$f \approx 0.316 Re^{-1/4}$$

Blasius correlation for Re < 10⁵

Moody Diagram

The famous chart relating friction factor f to Reynolds number Re and relative roughness ε/D. Essential for pipe flow calculations. Uses Colebrook-White equation for turbulent flow with roughness effects.

Boundary Layers

Near solid surfaces, viscous effects create a thin boundary layer where velocity varies from zero (at wall, no-slip) to the free-stream value. Outside this layer, the flow can often be treated as inviscid.

Flat Plate Boundary Layer (Blasius)

For laminar flow over a flat plate:

$$\delta = \frac{5x}{\sqrt{Re_x}} \quad \text{where} \quad Re_x = \frac{U_\infty x}{\nu}$$

Boundary layer thickness grows as √x. At higher Re_x, the layer is thinner.

Wall Shear Stress

$$\tau_w = 0.332 \rho U_\infty^2 Re_x^{-1/2}$$

Skin Friction Coefficient

$$C_f = \frac{0.664}{\sqrt{Re_x}}$$

Boundary Layer Thickness Definitions

δ (99%)

Distance where u = 0.99U_∞

δ* (Displacement)

Effect on outer flow: δ* = ∫(1-u/U)dy

θ (Momentum)

Momentum deficit: θ = ∫(u/U)(1-u/U)dy

Boundary Layer Separation

Adverse pressure gradients (dp/dx > 0) can cause flow reversal near the wall, leading to boundary layer separation. This creates wakes, increases drag dramatically, and can cause stall in airfoils.

Turbulence

Turbulent flows exhibit chaotic, irregular fluctuations with a wide range of length and time scales. Direct simulation (DNS) of all scales is often impractical, leading to modeling approaches.

Reynolds Decomposition

Split velocity into mean and fluctuating parts:

$$u = \bar{u} + u' \quad \text{where} \quad \bar{u'} = 0$$

Reynolds-Averaged Navier-Stokes (RANS)

Time-averaging introduces Reynolds stresses:

$$\rho\left(\frac{\partial \bar{u}_i}{\partial t} + \bar{u}_j\frac{\partial \bar{u}_i}{\partial x_j}\right) = -\frac{\partial \bar{p}}{\partial x_i} + \mu\nabla^2\bar{u}_i - \rho\frac{\partial \overline{u'_i u'_j}}{\partial x_j}$$

The Reynolds stress tensor -ρu'_iu'_j requires closure models (k-ε, k-ω, etc.)

Kolmogorov Scales

Smallest scales η = (ν³/ε)^1/4 where dissipation occurs. Energy cascades from large eddies to small ones.

Turbulent Viscosity

Eddy viscosity μ_t models enhanced mixing: τ_turb = μ_t ∂ū/∂y

Interactive Simulations

Poiseuille Flow: Velocity Profile and Flow Rate

Python

script.py46 lines

import numpy as np

# Hagen-Poiseuille flow in a circular pipe
# Fully developed, laminar, incompressible, steady
# u(r) = (dP/dx) * (R^2 - r^2) / (4*mu)

R = 0.025         # pipe radius (m)
L = 1.0           # pipe length (m)
mu = 1.002e-3     # dynamic viscosity of water at 20C (Pa.s)
rho = 998.0       # density (kg/m^3)
dP = 5000.0       # pressure drop (Pa)
dPdx = dP / L     # pressure gradient

print("Hagen-Poiseuille Flow in a Circular Pipe")
print("=" * 55)
print(f"Pipe radius: {R*1000:.1f} mm, Length: {L:.1f} m")
print(f"Viscosity: {mu:.3e} Pa.s, Pressure drop: {dP:.0f} Pa")
print()

# Velocity profile
r = np.linspace(0, R, 21)
u = (dPdx * (R**2 - r**2)) / (4 * mu)

print(f"{'r/R':>6} {'r (mm)':>8} {'u (m/s)':>10} {'u/u_max':>8}")
print("-" * 36)
u_max = u[0]
for ri, ui in zip(r, u):
    print(f"{ri/R:6.2f} {ri*1000:8.2f} {ui:10.4f} {ui/u_max:8.4f}")

# Flow rate (analytical)
Q = np.pi * R**4 * dPdx / (8 * mu)
u_avg = Q / (np.pi * R**2)
Re = rho * u_avg * (2*R) / mu

print(f"\nMax velocity (centerline): {u_max:.4f} m/s")
print(f"Average velocity: {u_avg:.4f} m/s")
print(f"u_avg / u_max = {u_avg/u_max:.4f} (theory: 0.5)")
print(f"Volume flow rate: {Q*1e6:.4f} mL/s")
print(f"Reynolds number: {Re:.1f}")
print(f"Flow regime: {'Laminar' if Re < 2300 else 'Turbulent'}")

# Wall shear stress
tau_w = dPdx * R / 2
print(f"Wall shear stress: {tau_w:.4f} Pa")
print(f"Darcy friction factor: {64/Re:.4f} (theory: 64/Re)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Reynolds Number and Boundary Layer Thickness

Fortran

program.f9057 lines

program boundary_layer
  implicit none
  ! Blasius boundary layer on a flat plate
  ! delta = 5.0 * x / sqrt(Re_x)
  ! delta* = 1.72 * x / sqrt(Re_x)  (displacement)
  ! theta = 0.664 * x / sqrt(Re_x)  (momentum)
  real(8) :: U_inf, nu, rho, L, x
  real(8) :: Re_x, Re_L, delta, delta_star, theta
  real(8) :: tau_w, Cf, Cd, F_drag
  real(8) :: W  ! plate width
  integer :: i

U_inf = 5.0d0     ! freestream velocity (m/s)
  rho = 1.225d0     ! air density (kg/m^3)
  nu = 1.5e-5       ! kinematic viscosity of air (m^2/s)
  L = 1.0d0         ! plate length (m)
  W = 0.5d0         ! plate width (m)

Re_L = U_inf * L / nu

write(*,'(A)') "Blasius Boundary Layer on a Flat Plate"
  write(*,'(A)') "======================================"
  write(*,'(A,F6.2,A)') "Freestream velocity: ", U_inf, " m/s"
  write(*,'(A,ES10.3)') "Reynolds number Re_L = ", Re_L
  write(*,*)
  write(*,'(A8,A12,A10,A10,A10,A10)') &
    "x(m)", "Re_x", "delta", "delta*", "theta", "Cf"
  write(*,'(A8,A12,A10,A10,A10,A10)') &
    "", "", "(mm)", "(mm)", "(mm)", ""
  write(*,'(A)') repeat("-", 62)

do i = 1, 10
    x = dble(i) * 0.1d0
    Re_x = U_inf * x / nu

if (Re_x > 0.0d0) then
      delta = 5.0d0 * x / sqrt(Re_x)
      delta_star = 1.72d0 * x / sqrt(Re_x)
      theta = 0.664d0 * x / sqrt(Re_x)
      Cf = 0.664d0 / sqrt(Re_x)

write(*,'(F8.2,ES12.3,F10.3,F10.3,F10.3,ES10.3)') &
        x, Re_x, delta*1000, delta_star*1000, &
        theta*1000, Cf
    end if
  end do

! Total drag on plate (both sides)
  Cd = 1.328d0 / sqrt(Re_L)
  F_drag = Cd * 0.5d0 * rho * U_inf**2 * L * W * 2.0d0

write(*,*)
  write(*,'(A,ES10.3)') "Total drag coefficient Cd = ", Cd
  write(*,'(A,F10.6,A)') "Total drag force (both sides): ", &
    F_drag, " N"
end program boundary_layer

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Detailed Derivations

Derivation: Plane Couette Flow

Consider flow between two infinite parallel plates separated by distance h. The bottom plate is stationary and the top plate moves with velocity U. Flow is steady, fully developed, incompressible, with no pressure gradient.

Simplify the x-momentum Navier-Stokes equation. Fully developed means ∂u/∂x = 0. Steady means ∂u/∂t = 0. No pressure gradient: ∂p/∂x = 0. No gravity in x. Only u(y) is nonzero:

$$0 = \mu\frac{d^2u}{dy^2}$$

Integrate twice:

$$\frac{du}{dy} = C_1 \implies u(y) = C_1\,y + C_2$$

Apply boundary conditions: u(0) = 0 (no-slip at bottom) and u(h) = U (no-slip at top):

$$C_2 = 0, \quad C_1 = \frac{U}{h}$$

Final velocity profile — linear:

$$\boxed{u(y) = U\,\frac{y}{h}}$$

The shear stress is uniform: τ = μU/h. This is the simplest viscous flow and models lubrication, bearings, and viscometry.

Derivation: Hagen-Poiseuille Flow in a Pipe

For steady, fully developed, axisymmetric flow in a pipe of radius R, the Navier-Stokes equation in cylindrical coordinates (r, θ, z) reduces to:

$$0 = -\frac{dp}{dz} + \mu\frac{1}{r}\frac{d}{dr}\!\left(r\frac{du}{dr}\right)$$

where u = u(r) is the axial velocity and dp/dz = −ΔP/L = const.

Rearrange and integrate once in r:

$$\frac{d}{dr}\!\left(r\frac{du}{dr}\right) = \frac{1}{\mu}\frac{dp}{dz}\,r \implies r\frac{du}{dr} = \frac{1}{2\mu}\frac{dp}{dz}\,r^2 + C_1$$

Boundedness at r = 0 requires C₁ = 0. Integrate again:

$$u(r) = \frac{1}{4\mu}\frac{dp}{dz}\,r^2 + C_2$$

No-slip at r = R: u(R) = 0. Substituting dp/dz = −ΔP/L:

$$\boxed{u(r) = \frac{\Delta P}{4\mu L}\left(R^2 - r^2\right)}$$

Parabolic velocity profile. Maximum at centreline: u_max = ΔP R²/(4μL).

Volume flow rate by integrating over the cross-section:

$$Q = \int_0^R u(r)\,2\pi r\,dr = \frac{\pi\Delta P}{8\mu L}\int_0^R (R^2 - r^2)\,r\,dr$$

Evaluating the integral:

$$\boxed{Q = \frac{\pi R^4 \Delta P}{8\mu L}}$$

Hagen-Poiseuille law. Flow rate is proportional to R⁴ — halving the radius reduces flow by a factor of 16.

Derivation: Stokes' Law for Creeping Flow Past a Sphere

For very slow flow (Re « 1, "creeping flow"), inertial terms in Navier-Stokes are negligible compared to viscous terms. The governing equations simplify to:

$$\nabla p = \mu\nabla^2\vec{v}, \quad \nabla\cdot\vec{v} = 0$$

For flow past a sphere of radius R moving at velocity U (or equivalently uniform flow U past a stationary sphere), assume the solution has axial symmetry. Use a stream function approach in spherical coordinates (r, θ):

$$\psi(r,\theta) = U\sin^2\!\theta\left(\frac{r^2}{2} - \frac{3Rr}{4} + \frac{R^3}{4r}\right)$$

The drag has two contributions. Integrating the pressure distribution over the sphere gives the form (pressure) drag. Integrating the shear stress gives the friction drag:

$$F_{\text{pressure}} = 2\pi\mu RU, \quad F_{\text{friction}} = 4\pi\mu RU$$

Total drag force (Stokes' law):

$$\boxed{F_D = 6\pi\mu R U}$$

Two-thirds of the drag is friction, one-third is pressure. The drag coefficient is C_D = 24/Re. Valid for Re < 1.

Derivation: Blasius Boundary Layer Solution

For steady laminar flow over a flat plate at zero incidence, the boundary layer equations (Prandtl, 1904) are:

$$u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = \nu\frac{\partial^2 u}{\partial y^2}, \quad \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0$$

Blasius introduced the similarity variable and stream function:

$$\eta = y\sqrt{\frac{U_\infty}{\nu x}}, \quad \psi = \sqrt{\nu x U_\infty}\,f(\eta)$$

Then u = U_∞f'(η) and the two PDEs collapse to a single ODE.

The Blasius equation (third-order nonlinear ODE):

$$\boxed{f''' + \frac{1}{2}f\,f'' = 0}$$

BCs: f(0) = 0, f'(0) = 0 (no-slip), f'(∞) = 1 (freestream).

Numerical solution shows u/U_∞ = 0.99 at η ≈ 5.0, giving the boundary layer thickness:

$$\boxed{\frac{\delta}{x} = \frac{5.0}{\sqrt{Re_x}}}$$

Displacement and momentum thicknesses:

$$\delta^* = \frac{1.72\,x}{\sqrt{Re_x}}, \quad \theta = \frac{0.664\,x}{\sqrt{Re_x}}$$

The wall shear stress: τ_w = 0.332 ρU_∞² Re_x^−1/2. Skin friction coefficient: C_f = 0.664/√Re_x.

Derivation: Critical Reynolds Number for Flow Transition

The Reynolds number represents the ratio of inertial to viscous forces. From the Navier-Stokes equations, non-dimensionalise using characteristic scales L, U, ρU² for length, velocity, and pressure:

$$\frac{\partial \vec{v}^*}{\partial t^*} + (\vec{v}^*\cdot\nabla^*)\vec{v}^* = -\nabla^* p^* + \frac{1}{Re}\nabla^{*2}\vec{v}^*$$

where Re = UL/ν. At low Re, the viscous term (1/Re)∇²v dominates; at high Re, the inertial convective term dominates.

Pipe flow (Osborne Reynolds, 1883): Reynolds' dye experiment showed laminar flow transitions to turbulent flow. The critical Reynolds number based on pipe diameter D and mean velocity V is:

$$\boxed{Re_{cr} = \frac{V D}{\nu} \approx 2300}$$

Below 2300: laminar. Above ~4000: fully turbulent. The range 2300–4000 is the transition region where intermittent turbulence (puffs, slugs) appears.

Flat plate boundary layer: Linear stability analysis of the Blasius profile shows that small disturbances (Tollmien-Schlichting waves) grow when:

$$\boxed{Re_{x,cr} = \frac{U_\infty x_{cr}}{\nu} \approx 5 \times 10^5}$$

This depends on freestream turbulence, surface roughness, and pressure gradients. In very quiet environments, transition can be delayed to Re_x > 10⁶.

Velocity Profiles: Couette, Poiseuille & Boundary Layer

Side-by-side comparison of the three fundamental viscous flow velocity profiles.

Advanced Simulations

Couette & Poiseuille Velocity Profiles

Python

script.py71 lines

import numpy as np

# === COUETTE FLOW ===
print("=" * 55)
print("=== PLANE COUETTE FLOW: u(y) = U * y/h ===")
print("=" * 55)

U = 2.0       # top plate velocity (m/s)
h = 0.01      # gap width (m)
mu = 0.1      # dynamic viscosity (Pa.s)

y_norm = np.linspace(0, 1, 11)  # y/h
u_couette = U * y_norm

print("  U = {:.1f} m/s, h = {:.0f} mm, mu = {:.3f} Pa.s".format(U, h*1000, mu))
print("  Shear stress tau = mu*U/h = {:.1f} Pa (uniform)".format(mu*U/h))
print()
print("  {:>6s}  {:>8s}  {:>8s}".format("y/h", "y (mm)", "u (m/s)"))
print("  " + "-" * 28)
for yn, uc in zip(y_norm, u_couette):
    print("  {:6.2f}  {:8.3f}  {:8.3f}".format(yn, yn*h*1000, uc))

# === POISEUILLE FLOW (Pipe) ===
print()
print("=" * 55)
print("=== HAGEN-POISEUILLE FLOW: u(r) = u_max*(1-r^2/R^2) ===")
print("=" * 55)

R = 0.025      # pipe radius (m)
dP = 10000.0   # pressure drop (Pa)
L = 1.0        # pipe length (m)
mu_w = 1.002e-3

u_max = dP * R**2 / (4 * mu_w * L)
r_norm = np.linspace(0, 1, 11)  # r/R
u_pois = u_max * (1 - r_norm**2)

Q = np.pi * R**4 * dP / (8 * mu_w * L)
V_avg = Q / (np.pi * R**2)
Re = 998 * V_avg * 2*R / mu_w

print("  R = {:.1f} mm, L = {:.1f} m, dP = {:.0f} Pa".format(R*1000, L, dP))
print("  u_max = {:.4f} m/s, V_avg = {:.4f} m/s (ratio = {:.2f})".format(
    u_max, V_avg, V_avg/u_max))
print("  Q = {:.4f} mL/s, Re = {:.1f}".format(Q*1e6, Re))
print()
print("  {:>6s}  {:>8s}  {:>10s}  {:>8s}".format("r/R", "r (mm)", "u (m/s)", "u/u_max"))
print("  " + "-" * 38)
for rn, up in zip(r_norm, u_pois):
    print("  {:6.2f}  {:8.2f}  {:10.4f}  {:8.4f}".format(rn, rn*R*1000, up, up/u_max))

# === COMBINED COUETTE-POISEUILLE (plane) ===
print()
print("=" * 55)
print("=== COUETTE-POISEUILLE FLOW (with pressure gradient) ===")
print("=" * 55)

# u(y) = U*y/h + (1/(2*mu))*(-dP/dx)*(hy - y^2)
dPdx_vals = [-5000, 0, 5000]  # Pa/m (favorable, zero, adverse)
print("  U = {:.1f} m/s, h = {:.0f} mm".format(U, h*1000))
print()
for dPdx in dPdx_vals:
    label = "Favorable" if dPdx < 0 else ("Zero" if dPdx == 0 else "Adverse")
    print("  dP/dx = {:.0f} Pa/m ({})".format(dPdx, label))
    print("  {:>6s}  {:>10s}".format("y/h", "u (m/s)"))
    for yn in [0.0, 0.25, 0.5, 0.75, 1.0]:
        y = yn * h
        u_cp = U*y/h + (1/(2*mu))*(-dPdx)*(h*y - y**2)
        print("  {:6.2f}  {:10.4f}".format(yn, u_cp))
    print()

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Blasius Boundary Layer Solution (Numerical)

Python

script.py87 lines

import numpy as np

# Solve the Blasius equation: f''' + 0.5*f*f'' = 0
# BCs: f(0)=0, f'(0)=0, f'(inf)=1
# Use shooting method with RK4

def blasius_rk4(f0_pp, eta_max=10.0, N=1000):
    """Integrate Blasius equation using RK4 with initial f''(0)=f0_pp."""
    deta = eta_max / N
    # State: [f, f', f'']
    y = np.array([0.0, 0.0, f0_pp])
    eta_arr = [0.0]
    f_arr = [y.copy()]

def deriv(y_loc):
        return np.array([y_loc[1], y_loc[2], -0.5*y_loc[0]*y_loc[2]])

for i in range(N):
        k1 = deriv(y)
        k2 = deriv(y + 0.5*deta*k1)
        k3 = deriv(y + 0.5*deta*k2)
        k4 = deriv(y + deta*k3)
        y = y + (deta/6)*(k1 + 2*k2 + 2*k3 + k4)
        eta_arr.append((i+1)*deta)
        f_arr.append(y.copy())

return np.array(eta_arr), np.array(f_arr)

# Shooting method: find f''(0) such that f'(inf) = 1
# Known answer: f''(0) = 0.33206
a, b = 0.3, 0.4
for _ in range(50):
    mid = (a + b) / 2
    eta, sol = blasius_rk4(mid)
    if sol[-1, 1] < 1.0:
        a = mid
    else:
        b = mid

f_pp_0 = (a + b) / 2
eta, sol = blasius_rk4(f_pp_0)

print("=== Blasius Boundary Layer Solution ===")
print("  f&apos;&apos;(0) = {:.5f} (exact: 0.33206)".format(f_pp_0))
print()
print("  {:>6s}  {:>10s}  {:>10s}  {:>10s}".format("eta", "f", "f&apos;=u/U", "f&apos;&apos;"))
print("  " + "-" * 42)
for i in range(0, len(eta), len(eta)//20):
    print("  {:6.2f}  {:10.5f}  {:10.5f}  {:10.5f}".format(
        eta[i], sol[i,0], sol[i,1], sol[i,2]))

# Find eta where f' = 0.99 (boundary layer edge)
for i in range(len(sol)):
    if sol[i,1] >= 0.99:
        print()
        print("  u/U = 0.99 at eta = {:.2f} => delta = {:.2f} * x/sqrt(Re_x)".format(
            eta[i], eta[i]))
        break

# Displacement and momentum thickness
deta = eta[1] - eta[0]
delta_star_coeff = np.trapezoid(1 - sol[:,1], eta)
theta_coeff = np.trapezoid(sol[:,1]*(1 - sol[:,1]), eta)
print("  Displacement thickness: delta*/x = {:.4f} / sqrt(Re_x)".format(delta_star_coeff))
print("  Momentum thickness:     theta/x  = {:.4f} / sqrt(Re_x)".format(theta_coeff))
print("  Shape factor H = delta*/theta = {:.4f}".format(delta_star_coeff/theta_coeff))

# Practical example
print()
print("=== Practical Example: Air over a 1m flat plate ===")
U_inf = 10.0
nu_air = 1.5e-5
rho_air = 1.225
x_vals = [0.1, 0.2, 0.5, 1.0]
print("  U = {:.0f} m/s, nu = {:.1e} m2/s".format(U_inf, nu_air))
print()
print("  {:>6s}  {:>10s}  {:>10s}  {:>10s}  {:>10s}".format(
    "x (m)", "Re_x", "delta(mm)", "Cf", "tau_w(Pa)"))
print("  " + "-" * 52)
for x in x_vals:
    Rex = U_inf * x / nu_air
    delta = 5.0 * x / np.sqrt(Rex)
    Cf = 0.664 / np.sqrt(Rex)
    tau_w = Cf * 0.5 * rho_air * U_inf**2
    print("  {:6.2f}  {:10.0f}  {:10.3f}  {:10.6f}  {:10.4f}".format(
        x, Rex, delta*1000, Cf, tau_w))

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Drag Coefficient vs Reynolds Number for a Sphere

Python

script.py88 lines

import numpy as np

# Drag coefficient correlations for a sphere across Re range
# Stokes: CD = 24/Re  (Re < 1)
# Oseen:  CD = 24/Re * (1 + 3Re/16)  (Re < 5)
# Standard drag curve + empirical fits

print("=== Sphere Drag Coefficient vs Reynolds Number ===")
print()

def cd_sphere(Re):
    """Empirical drag coefficient for a smooth sphere."""
    if Re < 0.01:
        return 24.0/Re
    elif Re < 1:
        return 24.0/Re * (1 + 3*Re/16)  # Oseen
    elif Re < 1000:
        return 24.0/Re * (1 + 0.15*Re**0.687)  # Schiller-Naumann
    elif Re < 2e5:
        return 0.44  # Newton regime
    else:
        return 0.1   # post-drag-crisis (turbulent BL)

Re_values = [0.001, 0.01, 0.1, 0.5, 1, 5, 10, 50, 100,
             500, 1000, 5000, 10000, 100000, 500000, 1000000]

print("  {:>12s}  {:>10s}  {:>12s}  {:>12s}".format(
    "Re", "CD", "CD*Re^2", "Regime"))
print("  " + "-" * 52)

for Re in Re_values:
    cd = cd_sphere(Re)
    if Re < 1:
        regime = "Stokes"
    elif Re < 1000:
        regime = "Intermediate"
    elif Re < 2e5:
        regime = "Newton"
    else:
        regime = "Drag crisis"
    print("  {:12.3f}  {:10.4f}  {:12.3f}  {:>12s}".format(
        Re, cd, cd*Re**2, regime))

# Terminal velocity of falling spheres
print()
print("=" * 55)
print("=== Terminal Velocity of Spheres Falling in Water ===")
print("=" * 55)

rho_f = 998.0
mu_f = 1.002e-3
g = 9.81

print("  Fluid: water (rho={:.0f}, mu={:.3e})".format(rho_f, mu_f))
print()

spheres = [
    ("Sand grain",    2650.0, 0.0005),
    ("Glass bead",    2500.0, 0.002),
    ("Steel ball",    7800.0, 0.005),
    ("Ping pong",       80.0, 0.020),
    ("Baseball",      700.0, 0.037),
]

print("  {:14s} {:>8s} {:>8s} {:>10s} {:>8s} {:>10s}".format(
    "Object", "rho", "D(mm)", "V_t(m/s)", "Re", "Regime"))
print("  " + "-" * 65)

for name, rho_s, D in spheres:
    R = D/2
    # Iterative: guess V, compute Re, get CD, compute new V
    V = 0.1
    for _ in range(200):
        Re = rho_f * V * D / mu_f
        cd = cd_sphere(max(Re, 0.001))
        V_new = np.sqrt(4*g*D*(rho_s - rho_f)/(3*cd*rho_f))
        if rho_s < rho_f:
            V_new = np.sqrt(4*g*D*(rho_f - rho_s)/(3*cd*rho_f))
        if abs(V_new - V) < 1e-8:
            break
        V = 0.5*(V + V_new)

Re_final = rho_f * V * D / mu_f
    regime = "Stokes" if Re_final < 1 else ("Interm." if Re_final < 1000 else "Newton")
    direction = "rises" if rho_s < rho_f else "falls"
    print("  {:14s} {:8.0f} {:8.1f} {:10.4f} {:8.0f} {:>10s}".format(
        name, rho_s, D*1000, V, int(Re_final), regime))

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Conservation Laws Video Lectures →