Computational Fluid Dynamics (CFD) Using MATLAB
Computational Fluid Dynamics (CFD) involves solving fluid flow equations using numerical methods. MATLAB provides a powerful environment for implementing CFD techniques, including discretization, visualization, and post-processing.
Introduction to CFD
CFD aims to solve the governing equations of fluid flow, such as the Navier-Stokes equations, which describe the motion of viscous fluid substances. These equations can be expressed as:
Continuity Equation:
\(\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0\)
Momentum Equation:
\(\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \mathbf{u}) = -\nabla p + \nabla \cdot (\mu \nabla \mathbf{u}) + \mathbf{f}\)
Example Problems in CFD using Matlab
Example 1: Basic CFD Problem: Solving 1D Heat Conduction
Consider a simple example of 1D steady-state heat conduction in a rod, governed by:
\(\frac{\mathrm{d}^2 T}{\mathrm{d} x^2} = 0\)
We discretize the domain using finite difference methods and solve it numerically:
% Parameters
L = 1; % Length of the rod
N = 10; % Number of nodes
dx = L / (N - 1); % Grid spacing
k = 100; % Thermal conductivity
T_left = 100; % Left boundary temperature
T_right = 200; % Right boundary temperature
% Initialize temperature array
T = zeros(N, 1);
T(1) = T_left;
T(end) = T_right;
% Construct coefficient matrix
A = diag(2 * ones(N-2, 1)) + diag(-1 * ones(N-3, 1), 1) + diag(-1 * ones(N-3, 1), -1);
b = zeros(N-2, 1);
b(1) = T_left;
b(end) = T_right;
% Solve for internal temperatures
T_internal = A \ b;
T(2:end-1) = T_internal;
disp(T); % Display the temperature distribution
The output represents the temperature distribution along the rod.
100.0
122.2
144.4
...
200.0
Example 2: 2D Steady-State Heat Conduction
The governing equation is:
\(\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} = 0\)
% Parameters
Nx = 10; Ny = 10; % Grid size
Lx = 1; Ly = 1; % Dimensions
dx = Lx / (Nx-1); dy = Ly / (Ny-1);
T = zeros(Nx, Ny); % Initialize temperature field
% Boundary conditions
T(:,1) = 100; % Left
T(:,end) = 200; % Right
T(1,:) = 150; % Top
T(end,:) = 50; % Bottom
% Iterative solution using Gauss-Seidel
for iter = 1:1000
for i = 2:Nx-1
for j = 2:Ny-1
T(i,j) = 0.25 * (T(i-1,j) + T(i+1,j) + T(i,j-1) + T(i,j+1));
end
end
end
disp(T); % Display temperature field
Sample output:
[100, ..., 200]
[150, ..., ...]
[...]
Advanced Example Problems in Computational Fluid Dynamics (CFD) Using MATLAB
Computational Fluid Dynamics (CFD) solves and analyzes fluid flow problems using numerical methods and algorithms. MATLAB provides versatile tools for CFD simulations, including partial differential equation solvers and visualization libraries.
Example 3: Solving the 2D Heat Equation
The 2D heat equation describes the distribution of temperature in a rectangular region over time:
∂T/∂t = α(∂²T/∂x² + ∂²T/∂y²)
\[ \frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right) \]
Here’s how to solve it using finite difference methods:
% Parameters
Lx = 1; Ly = 1; % Domain size
Nx = 20; Ny = 20; % Number of grid points
alpha = 0.01; % Thermal diffusivity
dt = 0.01; % Time step
T_end = 1; % Simulation time
% Discretization
dx = Lx / (Nx - 1);
dy = Ly / (Ny - 1);
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
% Initial condition
T = zeros(Nx, Ny);
T(round(Nx/2), round(Ny/2)) = 100; % Initial heat at the center
% Time-stepping loop
for t = 0:dt:T_end
T_new = T;
for i = 2:Nx-1
for j = 2:Ny-1
T_new(i, j) = T(i, j) + alpha * dt * ((T(i+1, j) - 2*T(i, j) + T(i-1, j)) / dx^2 + ...
(T(i, j+1) - 2*T(i, j) + T(i, j-1)) / dy^2);
end
end
T = T_new;
end
% Visualization
[X, Y] = meshgrid(x, y);
surf(X, Y, T');
xlabel('x'); ylabel('y'); zlabel('Temperature');
title('2D Heat Distribution');
This code simulates heat diffusion in a rectangular domain and visualizes the temperature distribution at the end.
Example 4: Simulating Laminar Flow in a 2D Channel
Solving the Navier-Stokes equations for incompressible flow in a 2D channel:
∂u/∂t + u∂u/∂x + v∂u/∂y = -∂p/∂x + ν(∂²u/∂x² + ∂²u/∂y²)
∂v/∂t + u∂v/∂x + v∂v/∂y = -∂p/∂y + ν(∂²v/∂x² + ∂²v/∂y²)
\[ \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = -\frac{\partial p}{\partial x} + \nu \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \]
\[ \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = -\frac{\partial p}{\partial y} + \nu \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right) \]
Here’s an example of a steady-state solution using a simple finite difference method:
% Parameters
Lx = 1; Ly = 0.2; % Channel size
Nx = 50; Ny = 10; % Grid points
Re = 100; % Reynolds number
dx = Lx / (Nx - 1);
dy = Ly / (Ny - 1);
% Initialization
u = zeros(Nx, Ny);
v = zeros(Nx, Ny);
p = zeros(Nx, Ny);
% Boundary conditions
u(:, end) = 1; % Top wall moving
u(:, 1) = 0; % Bottom wall
% Iterative solver
for iter = 1:5000
% Update velocity and pressure fields (simplified pseudo-code)
% Solve momentum equations and enforce continuity
end
% Visualization
[X, Y] = meshgrid(linspace(0, Lx, Nx), linspace(0, Ly, Ny));
quiver(X, Y, u', v');
xlabel('x'); ylabel('y');
title('Velocity Field in 2D Channel');
Example 5: Vortex Shedding Around a Cylinder
Simulating flow around a cylinder can demonstrate vortex shedding phenomena. Here’s a setup using the immersed boundary method:
% Parameters
D = 0.1; % Cylinder diameter
Re = 100; % Reynolds number
Lx = 1; Ly = 0.5; % Domain size
Nx = 100; Ny = 50; % Grid points
% Initialize flow variables
u = zeros(Nx, Ny); v = zeros(Nx, Ny);
p = zeros(Nx, Ny);
% Cylinder mask
[X, Y] = meshgrid(linspace(0, Lx, Nx), linspace(0, Ly, Ny));
cylinder = ((X - Lx/4).^2 + (Y - Ly/2).^2) < (D/2)^2;
% Time-stepping loop
for t = 1:1000
% Solve Navier-Stokes equations with immersed boundary
end
% Visualization
contourf(X, Y, sqrt(u.^2 + v.^2));
colorbar;
xlabel('x'); ylabel('y');
title('Vortex Shedding Around Cylinder');
Example 6: Solving the 1D Advection Equation
The 1D advection equation models the transport of a scalar quantity u(x,t) in a domain:
∂u/∂t + c∂u/∂x = 0
\[ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 \]
We solve this using the Lax-Wendroff scheme, a second-order accurate numerical method:
% Parameters
L = 1; % Domain length
Nx = 100; % Number of grid points
dx = L / Nx; % Grid spacing
c = 1; % Advection speed
dt = 0.01; % Time step
T = 2; % Total simulation time
% Discretization
x = linspace(0, L, Nx);
u = exp(-100 * (x - 0.5).^2); % Initial condition
% Lax-Wendroff Scheme
for t = 0:dt:T
u_new = u;
for i = 2:Nx-1
u_new(i) = u(i) - 0.5 * c * dt / dx * (u(i+1) - u(i-1)) + ...
0.5 * (c * dt / dx)^2 * (u(i+1) - 2*u(i) + u(i-1));
end
u = u_new;
end
% Plot results
plot(x, u);
xlabel('x'); ylabel('u');
title('1D Advection Equation Solution');
This simulates the transport of a Gaussian pulse in a 1D domain.
Example 7: Simulating Flow in a Lid-Driven Cavity
The lid-driven cavity problem involves solving the incompressible Navier-Stokes equations in a square domain with a moving top boundary.
% Parameters
L = 1; % Length of the cavity
N = 50; % Number of grid points
Re = 100; % Reynolds number
dx = L / (N-1);
dy = dx;
% Initialize velocity and pressure fields
u = zeros(N, N); v = zeros(N, N);
p = zeros(N, N);
% Boundary conditions
u(:, end) = 1; % Lid moves with unit velocity
% Iterative solver
for iter = 1:1000
% Solve momentum equations (simplified for brevity)
% Update pressure using the continuity equation
end
% Visualize velocity field
[X, Y] = meshgrid(linspace(0, L, N), linspace(0, L, N));
quiver(X, Y, u', v');
xlabel('x'); ylabel('y');
title('Lid-Driven Cavity Flow');
Example 8: Modeling Turbulent Channel Flow
Turbulent flows often require advanced techniques such as Reynolds-Averaged Navier-Stokes (RANS) equations or Large Eddy Simulation (LES). Here’s a simplified example of a RANS-based approach:
% Parameters
Lx = 2; Ly = 0.5; % Channel dimensions
Nx = 100; Ny = 50; % Grid points
Re_tau = 180; % Friction Reynolds number
% Initialize fields
u = zeros(Nx, Ny); % Streamwise velocity
k = 0.1 * ones(Nx, Ny); % Turbulence kinetic energy
epsilon = 0.01 * ones(Nx, Ny); % Turbulence dissipation rate
% Iterative solver
for iter = 1:5000
% Update k, epsilon, and velocity using RANS equations
end
% Visualization
[X, Y] = meshgrid(linspace(0, Lx, Nx), linspace(0, Ly, Ny));
contourf(X, Y, u');
colorbar;
xlabel('x'); ylabel('y');
title('Turbulent Channel Flow: Streamwise Velocity');
Example 9: Shock Wave Simulation Using the Euler Equations
The Euler equations describe compressible flows, often used for modeling shock waves:
∂ρ/∂t + ∂(ρu)/∂x = 0
∂(ρu)/∂t + ∂(ρu² + p)/∂x = 0
∂E/∂t + ∂((E+p)u)/∂x = 0
\[ \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x} \left( \rho u \right) = 0 \]
\[ \frac{\partial (\rho u)}{\partial t} + \frac{\partial (\rho u^2 + p)}{\partial x} = 0 \]
\[ \frac{\partial E}{\partial t} + \frac{\partial}{\partial x} \left( (E + p)u \right) = 0 \]
% Parameters
L = 1; Nx = 100; % Domain length and grid points
dx = L / Nx;
gamma = 1.4; % Specific heat ratio
% Initial conditions
rho = [ones(1, Nx/2), 0.125 * ones(1, Nx/2)];
u = zeros(1, Nx);
p = [ones(1, Nx/2), 0.1 * ones(1, Nx/2)];
% Conserved variables
E = p / (gamma - 1) + 0.5 * rho .* u.^2;
% Time-stepping loop
for t = 1:1000
% Update using numerical flux (e.g., Roe or HLL solver)
end
% Visualization
plot(linspace(0, L, Nx), rho);
xlabel('x'); ylabel('Density');
title('Shock Wave Simulation');
Example 10: Single-Stage Shock Tube Simulation
The single-stage shock tube problem involves a diaphragm separating two regions of gas at different pressures. At time t=0, the diaphragm is removed, and a shock wave propagates into the low-pressure region.
The governing equations are the Euler equations:
∂ρ/∂t + ∂(ρu)/∂x = 0
∂(ρu)/∂t + ∂(ρu² + p)/∂x = 0
∂E/∂t + ∂((E+p)u)/∂x = 0
% Parameters
gamma = 1.4; % Specific heat ratio
Nx = 200; % Number of grid points
L = 1.0; % Length of the tube
dx = L / Nx; % Spatial resolution
dt = 0.0005; % Time step
T = 0.2; % Total simulation time
% Initial conditions
x = linspace(0, L, Nx);
rho = [ones(1, Nx/2), 0.125 * ones(1, Nx/2)];
u = zeros(1, Nx);
p = [ones(1, Nx/2), 0.1 * ones(1, Nx/2)];
E = p / (gamma - 1) + 0.5 * rho .* u.^2; % Total energy
% Simulation loop
for t = 0:dt:T
% Compute fluxes (Roe or HLL solver)
% Update conserved variables: rho, rho*u, and E
end
% Visualization
subplot(3, 1, 1); plot(x, rho); title('Density');
subplot(3, 1, 2); plot(x, u); title('Velocity');
subplot(3, 1, 3); plot(x, p); title('Pressure');
This simulation displays the evolution of density, velocity, and pressure profiles over time.
Example 11: Double-Stage Shock Tube Simulation
The double-stage shock tube introduces an intermediate diaphragm, creating three regions with different initial pressures. At t=0, the diaphragms are removed sequentially or simultaneously, generating complex wave interactions.
The governing equations remain the same Euler equations as used in Example 10
% Parameters
gamma = 1.4; % Specific heat ratio
Nx = 300; % Number of grid points
L = 1.5; % Length of the tube
dx = L / Nx; % Spatial resolution
dt = 0.0005; % Time step
T = 0.3; % Total simulation time
% Initial conditions (Three regions)
x = linspace(0, L, Nx);
rho = [ones(1, Nx/3), 0.5 * ones(1, Nx/3), 0.125 * ones(1, Nx/3)];
u = zeros(1, Nx);
p = [1.0 * ones(1, Nx/3), 0.4 * ones(1, Nx/3), 0.1 * ones(1, Nx/3)];
E = p / (gamma - 1) + 0.5 * rho .* u.^2; % Total energy
% Simulation loop
for t = 0:dt:T
% Compute fluxes (Roe or HLL solver)
% Update conserved variables: rho, rho*u, and E
end
% Visualization
subplot(3, 1, 1); plot(x, rho); title('Density');
subplot(3, 1, 2); plot(x, u); title('Velocity');
subplot(3, 1, 3); plot(x, p); title('Pressure');
This simulation shows complex wave interactions due to the initial three-region setup, demonstrating rarefactions and shock reflections.
Useful MATLAB Functions for CFD
meshgrid generates a grid for 2D/3D computations.surf creates a 3D surface plot.contour generates 2D contour plots.quiver plots vector fields.laplacian computes Laplacians for discretized data.Practice Questions
Test Yourself
1. Implement a 1D unsteady heat conduction problem using an explicit scheme.
2. Solve a 2D Navier-Stokes problem for laminar flow in a channel.
3. Use MATLAB to visualize a 2D velocity field using the quiver function.
4. Simulate heat diffusion in a 3D domain and visualize the results using slice plots.
5. Implement and simulate the flow in a 2D square cavity with moving lid boundary conditions.
6. Solve the unsteady Burgers’ equation using an explicit finite difference method.
7. Implement and solve the 1D Burgers’ equation using explicit and implicit schemes. Compare results.
8. Simulate the evolution of a vortex pair in 2D using the streamfunction-vorticity formulation.
9. Model supersonic flow over a wedge and visualize shock wave patterns.
10. Modify the single-stage shock tube to include heat conduction effects. Analyze the changes in results.
11. Implement a TVD (Total Variation Diminishing) scheme for the double-stage shock tube simulation. Compare with standard results.
12. Simulate a three-dimensional shock tube problem and visualize the shock wave propagation in 3D.