IWTLP - Learn Programming by Doing

Computational Fluid Dynamics (CFD) uses numerical methods to solve the equations governing fluid flow. For aerospace engineers, CFD is indispensable for analyzing aerodynamic performance, thermal management, propulsion systems, and environmental control systems.

Governing Equations

All of CFD is built on three fundamental conservation laws expressed as the Navier-Stokes equations.

Conservation of Mass (Continuity)

\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0

For incompressible flow ( $\rho$ = constant):

\nabla \cdot \mathbf{u} = 0

Conservation of Momentum

\rho \frac{\partial \mathbf{u}}{\partial t} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} = -\nabla p + \mu \nabla^2 \mathbf{u} + \rho \mathbf{g}

This is Newton's second law applied to a fluid element: the left side is inertia (mass × acceleration), and the right side is the sum of pressure forces, viscous forces, and body forces.

Conservation of Energy

\rho c_p \frac{\partial T}{\partial t} + \rho c_p (\mathbf{u} \cdot \nabla) T = k \nabla^2 T + \Phi

Where $\Phi$ is the viscous dissipation function and $k$ is thermal conductivity.

Discretization Methods

To solve these continuous equations on a computer, we must discretize them.

Finite Difference Method (FDM)

Replaces derivatives with difference quotients on a structured grid:

Forward difference:

\frac{\partial u}{\partial x} \approx \frac{u_{i+1} - u_i}{\Delta x}

Central difference (2nd order):

\frac{\partial u}{\partial x} \approx \frac{u_{i+1} - u_{i-1}}{2\Delta x}

Second derivative:

\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2}

Finite Volume Method (FVM)

Integrates conservation equations over control volumes, naturally enforcing conservation:

\frac{\partial}{\partial t} \int_V \rho \phi \, dV + \oint_S \rho \phi \mathbf{u} \cdot d\mathbf{S} = \oint_S \Gamma \nabla \phi \cdot d\mathbf{S} + \int_V S_\phi \, dV

This integral form means: rate of change = convective flux + diffusive flux + source terms.

FVM is the most widely used method in aerospace CFD (OpenFOAM, ANSYS Fluent, Star-CCM+).

Finite Element Method (FEM)

Approximates the solution as a weighted sum of basis functions:

u(\mathbf{x}) \approx \sum_{j=1}^{N} u_j \phi_j(\mathbf{x})

Commonly used for structural analysis and multiphysics problems. COMSOL and FEniCS use FEM.

Mesh Generation

Structured Meshes

Regular grid with implicit connectivity (i, j, k indexing):

Advantages: Simple data structure, efficient solvers, high accuracy Disadvantages: Difficult for complex geometries, mesh quality can degrade

Unstructured Meshes

Triangles (2D) or tetrahedra (3D) with explicit connectivity:

Advantages: Handles complex geometries, easy local refinement Disadvantages: More memory, slower solver, more complex code

Mesh Quality Metrics

Metric	Ideal	Acceptable
Aspect ratio	1.0	< 5.0
Skewness	0.0	< 0.85
Orthogonality	1.0	> 0.15
Expansion ratio	1.0	< 1.3

Wall Resolution

For boundary layer flows, the first cell height must satisfy:

y^+ = \frac{u_\tau y}{\nu} \approx 1 \quad \text{(for resolved BL)}

Where $u_\tau = \sqrt{\tau_w / \rho}$ is the friction velocity. This typically requires very thin cells near walls ( $y \sim 10^{-5}$ m for typical Reynolds numbers).

Time Integration

Explicit Methods

Forward Euler:

u^{n+1} = u^n + \Delta t \cdot f(u^n, t^n)

Simple but conditionally stable. Must satisfy the CFL condition:

CFL = \frac{u \Delta t}{\Delta x} \leq 1

Runge-Kutta (4th order): Higher accuracy with four intermediate evaluations per time step. Still explicit, still CFL-limited.

Implicit Methods

Backward Euler:

u^{n+1} = u^n + \Delta t \cdot f(u^{n+1}, t^{n+1})

Unconditionally stable (no CFL restriction), but requires solving a system of equations each step. Better for steady-state problems and stiff systems.

Turbulence Modeling

Most aerospace flows are turbulent. Direct Numerical Simulation (DNS) resolves all scales but is computationally prohibitive for practical Reynolds numbers.

Reynolds-Averaged Navier-Stokes (RANS)

Decomposes velocity into mean and fluctuating components:

u = \bar{u} + u'

The Reynolds stress tensor $\overline{u_i' u_j'}$ requires modeling. Common models:

Spalart-Allmaras: One-equation model, good for attached aerospace flows
k-ε: Two-equation model, robust but poor near walls
k-ω SST: Two-equation model, excellent for aerospace boundary layers (Menter's Shear Stress Transport)

Large Eddy Simulation (LES)

Resolves large-scale eddies, models small scales with subgrid-scale (SGS) models:

\bar{\tau}_{ij} - \frac{1}{3}\bar{\tau}_{kk}\delta_{ij} = -2\nu_{SGS}\bar{S}_{ij}

More accurate than RANS but 10-100x more expensive. Used for unsteady flows, acoustics, and combustion.

Boundary Conditions

Common Aerospace BCs

Type	Application	Specification
Inlet	Freestream, engine inlet	Velocity or total pressure + direction
Outlet	Downstream boundary	Static pressure or zero-gradient
Wall (no-slip)	Aircraft surface	$\mathbf{u} = 0$ at wall
Symmetry	Half-model simulation	Zero normal velocity, zero normal gradients
Far-field	External aerodynamics	Characteristic-based (non-reflecting)
Periodic	Turbine blade passages	Matched flow at periodic boundaries

Verification and Validation

Grid Convergence Study

Run simulations on progressively refined meshes and check if the solution converges:

GCI_{fine} = \frac{F_s |e_{rel}|}{r^p - 1}

Where $F_s$ is a safety factor (typically 1.25), $e_{rel}$ is the relative change between grids, $r$ is the refinement ratio, and $p$ is the observed order of accuracy.

Richardson Extrapolation

Estimate the exact solution from two mesh levels:

f_{exact} \approx f_1 + \frac{f_1 - f_2}{r^p - 1}

Where $f_1$ is the fine-grid solution and $f_2$ is the coarse-grid solution.

Your Challenge: 1D Heat Equation Solver

Implement a 1D heat equation solver using finite differences:

\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}

# 1D Heat Equation Solver (Explicit Forward Euler)
import math

# Domain parameters
L = 1.0          # Length of rod (m)
nx = 50          # Number of spatial points
dx = L / (nx - 1)
alpha = 0.01     # Thermal diffusivity (m^2/s)

# Time parameters (CFL-limited)
dt = 0.4 * dx**2 / alpha  # Satisfy stability criterion
nt = 200         # Number of time steps

# Initial condition: sine wave
T = [math.sin(math.pi * i * dx / L) for i in range(nx)]
T[0] = 0.0       # Left boundary (fixed temperature)
T[-1] = 0.0      # Right boundary (fixed temperature)

# Time integration
for n in range(nt):
    T_new = T[:]
    for i in range(1, nx - 1):
        T_new[i] = T[i] + alpha * dt / dx**2 * (T[i+1] - 2*T[i] + T[i-1])
    T = T_new

# Analytical solution for comparison
T_exact = [math.exp(-alpha * (math.pi/L)**2 * nt*dt) * math.sin(math.pi * i*dx/L) for i in range(nx)]

# Calculate error
max_error = max(abs(T[i] - T_exact[i]) for i in range(nx))
print(f"Grid spacing: dx = {dx:.4f}")
print(f"Time step: dt = {dt:.6f}")
print(f"CFL number: {alpha * dt / dx**2:.3f}")
print(f"Final time: {nt * dt:.4f} s")
print(f"Max error vs analytical: {max_error:.6f}")
print(f"Temperature at center: {T[nx//2]:.6f} (exact: {T_exact[nx//2]:.6f})")

Extend this to 2D using the Alternating Direction Implicit (ADI) method. What happens if you violate the CFL condition? Try increasing dt and observe the instability.

Computational Fluid Dynamics (Python)