Conservation Laws

Mass, momentum, and energy - the fundamental equations of fluid motion

Introduction

The motion of fluids is governed by conservation laws - mass, momentum, and energy cannot be created or destroyed, only transported and converted. These fundamental principles, combined with constitutive relations, give us the complete equations of fluid dynamics.

We present both integral forms (useful for control volume analysis) and differential forms (useful for detailed flow field solutions).

Conservation of Mass (Continuity)

Integral Form

The rate of change of mass in a control volume equals the net mass flux through the surface:

$$\frac{\partial}{\partial t}\int_{CV} \rho \, dV + \oint_{CS} \rho (\vec{v} \cdot \hat{n}) \, dA = 0$$

For steady flow: ∮ρ(v·n)dA = 0 → Mass in = Mass out

Differential Form

Using the divergence theorem:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{v}) = 0$$

Or equivalently, using the material derivative:

$$\frac{D\rho}{Dt} + \rho (\nabla \cdot \vec{v}) = 0$$

Incompressible Flow

When Dρ/Dt = 0 (density of a fluid particle doesn't change), continuity simplifies to:

$$\nabla \cdot \vec{v} = 0$$

Valid for liquids and low-speed (Ma < 0.3) gas flows.

Conservation of Momentum

Integral Form (Control Volume)

Newton's second law applied to a control volume:

$$\sum \vec{F} = \frac{\partial}{\partial t}\int_{CV} \vec{v}\rho \, dV + \oint_{CS} \vec{v}\rho (\vec{v} \cdot \hat{n}) \, dA$$

Forces include pressure, gravity, viscous stresses, and external forces.

Cauchy Momentum Equation

General differential form (valid for any continuous medium):

$$\rho \frac{D\vec{v}}{Dt} = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \rho \vec{g}$$

where τ is the deviatoric stress tensor (viscous stresses for Newtonian fluids).

Navier-Stokes Equations

For incompressible Newtonian fluids with constant viscosity:

$$\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}$$

∂v/∂t

Unsteady

(v·∇)v

Convection

μ∇²v

Diffusion

-∇p + ρg

Pressure + Gravity

Euler Equations (Inviscid)

When viscous effects are negligible (high Re, away from boundaries):

$$\rho \frac{D\vec{v}}{Dt} = -\nabla p + \rho \vec{g}$$

Conservation of Energy

First Law for a Control Volume

Energy balance including work and heat transfer:

$$\dot{Q} - \dot{W}_s = \frac{\partial}{\partial t}\int_{CV} e\rho \, dV + \oint_{CS} \left(e + \frac{p}{\rho}\right)\rho (\vec{v} \cdot \hat{n}) \, dA$$

where e = u + v²/2 + gz is specific energy (internal + kinetic + potential)

Differential Energy Equation

For a compressible fluid:

$$\rho \frac{De}{Dt} = -\nabla \cdot \vec{q} - p(\nabla \cdot \vec{v}) + \boldsymbol{\tau}:\nabla\vec{v}$$

Terms: heat conduction, pressure work, viscous dissipation

For Incompressible Flow

Often temperature can be treated separately from mechanics (weak coupling), and the energy equation simplifies to tracking temperature evolution via heat conduction and advection.

Bernoulli Equation

A special case of energy/momentum conservation for steady, inviscid, incompressible flow along a streamline:

$$p + \frac{1}{2}\rho v^2 + \rho g z = \text{constant along streamline}$$

Static pressure

½ρv²

Dynamic pressure

ρgz

Hydrostatic pressure

Extended Forms

Unsteady Bernoulli

$$\frac{p}{\rho} + \frac{v^2}{2} + gz + \int \frac{\partial v}{\partial t} ds = \text{constant}$$

With Head Loss (Engineering)

$$\frac{p_1}{\rho g} + \frac{v_1^2}{2g} + z_1 = \frac{p_2}{\rho g} + \frac{v_2^2}{2g} + z_2 + h_L$$

h_L accounts for viscous losses (friction, minor losses)

Applications

Pitot Tube

Velocity from Bernoulli: v = √(2(p₀-p)/ρ) where p₀ is stagnation pressure

Venturi Meter

Flow rate from pressure difference using continuity + Bernoulli

Siphon

Flow driven by elevation difference - must stay above vapor pressure

Jet Propulsion

Momentum flux from nozzle provides thrust: F = ṁ(v_exit - v_inlet)

Interactive Simulations

Bernoulli Equation: Flow Through a Converging Nozzle

Python

script.py44 lines

import numpy as np

# Bernoulli equation along a streamline (inviscid, incompressible):
# P1 + 0.5*rho*v1^2 + rho*g*z1 = P2 + 0.5*rho*v2^2 + rho*g*z2

rho = 1000.0    # water density (kg/m^3)
g = 9.81        # gravitational acceleration (m/s^2)

# Converging nozzle (horizontal, z1 = z2)
D_inlet = 0.10   # inlet diameter (m)
P_inlet = 200000  # inlet pressure (Pa, gauge)
v_inlet = 2.0    # inlet velocity (m/s)
A_inlet = np.pi * (D_inlet/2)**2

print("Bernoulli Equation: Converging Nozzle")
print("=" * 55)
print(f"Inlet: D = {D_inlet*100:.1f} cm, v = {v_inlet:.1f} m/s, P = {P_inlet/1000:.1f} kPa")
print()

# Calculate flow for different outlet diameters
D_outlets = np.array([0.08, 0.06, 0.04, 0.03, 0.02])

print(f"{'D_out (cm)':>10} {'v_out (m/s)':>12} {'P_out (kPa)':>12} {'Q (L/s)':>10}")
print("-" * 48)

Q = A_inlet * v_inlet  # volume flow rate (continuity)

for D_out in D_outlets:
    A_out = np.pi * (D_out/2)**2
    v_out = Q / A_out              # from continuity
    # Bernoulli: P_out = P_inlet + 0.5*rho*(v_inlet^2 - v_out^2)
    P_out = P_inlet + 0.5 * rho * (v_inlet**2 - v_out**2)
    print(f"{D_out*100:10.1f} {v_out:12.2f} {P_out/1000:12.2f} {Q*1000:10.3f}")

print(f"\nVolume flow rate Q = {Q*1000:.3f} L/s (conserved)")
print(f"\nNote: As the nozzle narrows, velocity increases")
print(f"and pressure decreases (Venturi effect).")

# Find cavitation limit
P_vapor = 2340  # water vapor pressure at 20C
v_cav = np.sqrt(2*(P_inlet - P_vapor)/rho + v_inlet**2)
D_cav = np.sqrt(4*Q/(np.pi*v_cav))
print(f"\nCavitation onset at D = {D_cav*100:.2f} cm (v = {v_cav:.1f} m/s)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Mass and Momentum Conservation in a Pipe Junction

Fortran

program.f9061 lines

program pipe_junction
  implicit none
  ! Conservation laws at a T-junction:
  ! Mass: rho*A1*v1 = rho*A2*v2 + rho*A3*v3
  ! Momentum balance in x-direction
  real(8) :: pi, rho, D1, D2, D3
  real(8) :: A1, A2, A3, v1, v2, v3
  real(8) :: Q1, Q2, Q3, P1, P2, P3
  real(8) :: mdot1, mdot2, mdot3
  real(8) :: F_x, F_y
  real(8) :: split_ratio
  integer :: i

pi = 4.0d0 * atan(1.0d0)
  rho = 1000.0d0   ! water

D1 = 0.10d0  ! inlet pipe diameter (m)
  D2 = 0.08d0  ! straight outlet diameter (m)
  D3 = 0.06d0  ! branch outlet diameter (m)

A1 = pi * (D1/2.0d0)**2
  A2 = pi * (D2/2.0d0)**2
  A3 = pi * (D3/2.0d0)**2

v1 = 3.0d0   ! inlet velocity (m/s)
  Q1 = A1 * v1
  P1 = 150000.0d0  ! inlet pressure (Pa gauge)

write(*,'(A)') "Pipe T-Junction: Conservation Analysis"
  write(*,'(A)') "======================================="
  write(*,'(A,F6.1,A)') "Inlet: D=", D1*100, " cm"
  write(*,'(A,F6.1,A)') "Straight: D=", D2*100, " cm"
  write(*,'(A,F6.1,A)') "Branch: D=", D3*100, " cm"
  write(*,'(A,F6.2,A)') "Inlet velocity: ", v1, " m/s"
  write(*,*)
  write(*,'(A8,A10,A10,A10,A10,A10)') &
    "Split%", "Q2(L/s)", "Q3(L/s)", "v2(m/s)", "v3(m/s)", "mdot err"
  write(*,'(A)') repeat("-", 60)

do i = 1, 9
    split_ratio = dble(i) * 0.1d0
    Q2 = split_ratio * Q1
    Q3 = Q1 - Q2    ! mass conservation

v2 = Q2 / A2
    v3 = Q3 / A3

mdot1 = rho * Q1
    mdot2 = rho * Q2
    mdot3 = rho * Q3

write(*,'(F7.0,A,F10.3,F10.3,F10.2,F10.2,ES10.2)') &
      split_ratio*100, "%", Q2*1000, Q3*1000, &
      v2, v3, mdot1-mdot2-mdot3
  end do

write(*,*)
  write(*,'(A,F8.3,A)') "Total flow rate: ", Q1*1000, " L/s"
  write(*,'(A,F8.2,A)') "Total mass flow: ", rho*Q1, " kg/s"
end program pipe_junction

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Detailed Derivations

Derivation: Continuity Equation from Mass Conservation

Consider a fixed, infinitesimal control volume dx dy dz in a flow field. Mass is conserved: the rate of mass increase inside equals the net mass flux inward through the surfaces.

Rate of mass increase inside the element:

$$\frac{\partial \rho}{\partial t}\,dx\,dy\,dz$$

Net mass flux through x-faces: mass in at x minus mass out at x+dx:

$$[\rho u]_x\,dy\,dz - [\rho u]_{x+dx}\,dy\,dz = -\frac{\partial(\rho u)}{\partial x}\,dx\,dy\,dz$$

Similarly for y-faces and z-faces.

Set rate of accumulation equal to net influx, cancel the volume element:

$$\boxed{\frac{\partial \rho}{\partial t} + \frac{\partial(\rho u)}{\partial x} + \frac{\partial(\rho v)}{\partial y} + \frac{\partial(\rho w)}{\partial z} = 0}$$

In vector notation:

$$\frac{\partial \rho}{\partial t} + \nabla\cdot(\rho\vec{v}) = 0$$

For incompressible flow (ρ = const): ∇·v = 0.

Derivation: Navier-Stokes Equations from Newton's Second Law

Apply Newton's second law to an infinitesimal fluid element of mass dm = ρ dx dy dz:

$$dm\,\frac{D\vec{v}}{Dt} = d\vec{F}_{\text{body}} + d\vec{F}_{\text{surface}}$$

Body forces (gravity): dF_body = ρg dV.

Surface forces (pressure + viscous stress): the net force in the x-direction from the stress tensor σ_ij:

$$dF_x = \left(\frac{\partial \sigma_{xx}}{\partial x} + \frac{\partial \sigma_{yx}}{\partial y} + \frac{\partial \sigma_{zx}}{\partial z}\right)dx\,dy\,dz$$

Decompose the stress tensor into pressure and deviatoric (viscous) parts: σ_ij = −pδ_ij + τ_ij. For a Newtonian fluid:

$$\tau_{ij} = \mu\left(\frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}\right) + \lambda(\nabla\cdot\vec{v})\delta_{ij}$$

where λ = −2μ/3 (Stokes hypothesis).

Combining and dividing by dV, for incompressible flow (∇·v = 0):

$$\boxed{\rho\frac{D\vec{v}}{Dt} = -\nabla p + \mu\nabla^2\vec{v} + \rho\vec{g}}$$

The x-component written fully:

$$\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$$

Three coupled, nonlinear, second-order PDEs plus the continuity constraint — together they are four equations for four unknowns (u, v, w, p).

Derivation: Bernoulli Equation Along a Streamline

Start with the Euler equation (inviscid, μ = 0) for steady flow (∂/∂t = 0):

$$(\vec{v}\cdot\nabla)\vec{v} = -\frac{1}{\rho}\nabla p + \vec{g}$$

Project along a streamline direction ds (where ds is tangent to v). Let v = |v| and note that (v·∇)v · ds = v dv:

$$v\,dv = -\frac{1}{\rho}\,dp - g\,dz$$

For incompressible flow (ρ = const), integrate along the streamline:

$$\int v\,dv = -\frac{1}{\rho}\int dp - g\int dz$$

Carrying out the integration:

$$\boxed{p + \frac{1}{2}\rho v^2 + \rho g z = \text{constant along a streamline}}$$

Assumptions: (i) Steady flow, (ii) Inviscid (frictionless), (iii) Incompressible (ρ = const), (iv) Along a single streamline, (v) No shaft work or heat transfer. Violations of any assumption require corrections (head loss, unsteady term, etc.).

Derivation: Energy Equation from the First Law of Thermodynamics

The first law for a fluid element (specific internal energy e, specific enthalpy h = e + p/ρ):

$$\rho\frac{De}{Dt} = \text{(heat added)} + \text{(work done on element)}$$

Heat flux by conduction (Fourier's law): q = −k∇T. Net heat added to the element per unit volume:

$$-\nabla\cdot\vec{q} = \nabla\cdot(k\nabla T)$$

Work by surface stresses on the element: −p(∇·v) (compression work) + Φ (viscous dissipation):

$$\Phi = \tau_{ij}\frac{\partial v_i}{\partial x_j} = \mu\left(\frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}\right)\frac{\partial v_i}{\partial x_j} \geq 0$$

Viscous dissipation is always non-negative — it converts kinetic energy to heat.

Combining all terms:

$$\boxed{\rho\frac{De}{Dt} = -p(\nabla\cdot\vec{v}) + \Phi + \nabla\cdot(k\nabla T)}$$

For incompressible flow, ∇·v = 0 and the compression work term vanishes.

Derivation: Darcy-Weisbach Equation from Momentum Integral

Consider fully developed, steady flow in a horizontal pipe of diameter D and length L. The flow is driven by a pressure difference ΔP = P₁ − P₂.

Apply the momentum equation to a cylindrical control volume. For fully developed flow, momentum flux in equals momentum flux out. The balance is between pressure forces and wall shear:

$$\Delta P \cdot \frac{\pi D^2}{4} = \tau_w \cdot \pi D L$$

Solve for the pressure drop:

$$\Delta P = \frac{4\tau_w L}{D}$$

Define the Darcy friction factor f via the non-dimensionalisation τ_w = (f/8)ρV²:

$$\boxed{\Delta P = f\,\frac{L}{D}\,\frac{\rho V^2}{2}}$$

For laminar flow, comparison with the Hagen-Poiseuille solution gives f = 64/Re. For turbulent flow, f depends on Re and relative roughness ε/D — determined from the Colebrook equation or the Moody diagram.

Control Volume Diagram

Mass flux, pressure forces, and body forces acting on a fixed control volume in a flow field.

Advanced Simulations

Bernoulli Applications: Venturi Meter & Pitot Tube

Python

script.py51 lines

import numpy as np

rho = 1000.0   # water density (kg/m3)
g = 9.81

print("=" * 60)
print("=== VENTURI METER: Flow Rate from Pressure Difference ===")
print("=" * 60)

D1 = 0.10        # pipe diameter (m)
D2 = 0.04        # throat diameter (m)
A1 = np.pi/4 * D1**2
A2 = np.pi/4 * D2**2
Cd = 0.98         # discharge coefficient

dP_values = [5000, 10000, 20000, 50000, 100000]  # Pa

print("  D_pipe = {:.0f} mm, D_throat = {:.0f} mm, Cd = {:.2f}".format(D1*1000, D2*1000, Cd))
print()
print("  {:>12s}  {:>10s}  {:>10s}  {:>10s}".format("dP (kPa)", "V1 (m/s)", "V2 (m/s)", "Q (L/s)"))
print("  " + "-" * 48)

for dP in dP_values:
    # Bernoulli + continuity: V2 = sqrt(2*dP / (rho*(1 - (A2/A1)^2)))
    V2 = Cd * np.sqrt(2*dP / (rho * (1 - (A2/A1)**2)))
    V1 = V2 * A2/A1
    Q = A2 * V2
    print("  {:12.1f}  {:10.3f}  {:10.3f}  {:10.3f}".format(dP/1000, V1, V2, Q*1000))

print()
print("=" * 60)
print("=== PITOT TUBE: Velocity from Stagnation Pressure ===")
print("=" * 60)

# Pitot tube: P_stag = P_static + 0.5*rho*V^2
# => V = sqrt(2*(P_stag - P_static)/rho)

rho_air = 1.225   # air density at sea level (kg/m3)
print("  Fluid: Air at sea level (rho = {:.3f} kg/m3)".format(rho_air))
print()
print("  {:>12s}  {:>12s}  {:>12s}".format("dP (Pa)", "V (m/s)", "V (km/h)"))
print("  " + "-" * 40)

for dP in [10, 50, 100, 500, 1000, 2000, 5000]:
    V = np.sqrt(2*dP/rho_air)
    print("  {:12.0f}  {:12.2f}  {:12.1f}".format(dP, V, V*3.6))

print()
print("  Note: Pitot tubes measure total (stagnation) pressure.")
print("  Combined with static ports, they give airspeed in aircraft.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Pipe Flow Pressure Drop & Moody Diagram Data

Python

script.py76 lines

import numpy as np

# Darcy-Weisbach: dP = f * (L/D) * (rho*V^2/2)
# Laminar: f = 64/Re
# Turbulent smooth: Blasius f = 0.316 * Re^(-1/4) for Re < 1e5
# Turbulent rough: Colebrook-White (iterative)

rho = 998.0       # water at 20C
mu = 1.002e-3     # dynamic viscosity (Pa.s)
nu = mu/rho
D = 0.05          # pipe diameter (m)
L = 100.0         # pipe length (m)
g = 9.81

print("=== Pipe Flow Pressure Drop Analysis ===")
print("  D = {:.0f} mm, L = {:.0f} m, fluid = water at 20C".format(D*1000, L))
print()

velocities = [0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0]
print("  {:>8s}  {:>10s}  {:>8s}  {:>12s}  {:>12s}".format(
    "V (m/s)", "Re", "f", "dP (kPa)", "Regime"))
print("  " + "-" * 56)

for V in velocities:
    Re = rho * V * D / mu
    if Re < 2300:
        f = 64.0 / Re
        regime = "Laminar"
    else:
        f = 0.316 * Re**(-0.25)  # Blasius (smooth pipe)
        regime = "Turbulent"
    dP = f * (L/D) * (rho * V**2 / 2)
    print("  {:8.3f}  {:10.0f}  {:8.5f}  {:12.3f}  {:>12s}".format(V, Re, f, dP/1000, regime))

# Moody diagram data: f vs Re for different roughnesses
print()
print("=" * 60)
print("=== Moody Diagram: f vs Re for Various Roughnesses ===")
print("=" * 60)
print()

def colebrook(Re, eD):
    """Solve Colebrook-White equation iteratively."""
    if Re < 2300:
        return 64.0/Re
    f = 0.02  # initial guess
    for _ in range(50):
        rhs = -2.0 * np.log10(eD/3.7 + 2.51/(Re*np.sqrt(f)))
        f_new = 1.0 / rhs**2
        if abs(f_new - f) < 1e-8:
            break
        f = f_new
    return f

roughness_ratios = [0, 1e-5, 1e-4, 1e-3, 5e-3]
Re_values = [500, 1000, 2000, 4000, 10000, 50000, 100000, 500000, 1000000]

header = "  {:>10s}".format("Re")
for eD in roughness_ratios:
    header += "  e/D={:.0e}".format(eD) if eD > 0 else "  {:>10s}".format("Smooth")
print(header)
print("  " + "-" * (12 + 12*len(roughness_ratios)))

for Re in Re_values:
    row = "  {:10.0f}".format(Re)
    for eD in roughness_ratios:
        f = colebrook(Re, eD)
        row += "  {:10.5f}".format(f)
    print(row)

print()
print("  Key observations:")
print("  - Laminar regime (Re<2300): f = 64/Re, independent of roughness")
print("  - Turbulent regime: f depends on both Re and e/D")
print("  - At high Re, f becomes independent of Re (fully rough regime)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Fluid Kinematics Viscous Flow →