Simulating Flow Over a Cylinder Using MATLAB and the Lattice Boltzmann Method

Introduction

Computational Fluid Dynamics (CFD) is a powerful tool used to analyze fluid flow problems. In this tutorial, we’ll simulate 2D flow past a cylinder using MATLAB and the Lattice Boltzmann Method (LBM) to compute the drag coefficient. We’ll also visualize velocity fields and streamline patterns.

Problem Setup

We’ll analyze the flow past a cylinder with the following parameters:

% Problem Parameters
D = 1.0;            % Cylinder diameter (m)
U = 1.0;            % Free-stream velocity (m/s)
rho = 1.0;          % Fluid density (kg/m^3)
mu = 0.01;          % Dynamic viscosity (Pa.s)
Re = rho * U * D / mu; % Reynolds number

% Domain size
Lx = 20;  % Length of domain in x-direction
Ly = 10;  % Height of domain in y-direction

fprintf('Reynolds number: %.2f\n', Re);

The Reynolds number provides insight into flow characteristics, helping to determine if the flow is laminar or turbulent.

Mesh Generation

We’ll create a rectangular computational domain and place the cylinder at the center:

% Grid resolution
Nx = 200; % Grid points in x-direction
Ny = 100; % Grid points in y-direction

x = linspace(-Lx/2, Lx/2, Nx);
y = linspace(-Ly/2, Ly/2, Ny);
[X, Y] = meshgrid(x, y);

% Define cylinder boundary condition
[theta, r] = cart2pol(X, Y);
cylinder_mask = r < D/2;   % Creates the surface mask

% Visualize the mesh
figure;
contourf(X, Y, double(cylinder_mask), [0.5, 0.5], 'k');
axis equal;
title('Mesh and Cylinder Location');
xlabel('X');
ylabel('Y');
grid on;

Solving the Navier-Stokes Equations Using Lattice Boltzmann Method

We’ll now solve the Navier-Stokes equations using the LBM approach with finite difference approximations.

% Initialize velocity and pressure fields
u = zeros(Ny, Nx);
v = zeros(Ny, Nx);
p = zeros(Ny, Nx);

dx = Lx / (Nx - 1);
dy = Ly / (Ny - 1);
dt = 0.01;  % Time step

maxIter = 1000;
for iter = 1:maxIter
    % Compute velocity with central difference scheme
    u_new = u - dt * ((u(:, [2:end, end]) - u(:, [1, 1:end-1])) / (2 * dx) ...
                      + (v([2:end, end], :) - v([1, 1:end-1], :)) / (2 * dy)) ...
            + dt * mu * ((u(:, [2:end, end]) - 2*u + u(:, [1, 1:end-1])) / dx^2 ...
                      + (u([2:end, end], :) - 2*u + u([1, 1:end-1], :)) / dy^2);

    v_new = v - dt * ((u(:, [2:end, end]) - u(:, [1, 1:end-1])) / (2 * dx) ...
                      + (v([2:end, end], :) - v([1, 1:end-1], :)) / (2 * dy)) ...
            + dt * mu * ((v(:, [2:end, end]) - 2*v + v(:, [1, 1:end-1])) / dx^2 ...
                      + (v([2:end, end], :) - 2*v + v([1, 1:end-1], :)) / dy^2);

    % Apply boundary conditions
    u_new(:, 1) = U;  
    u_new(:, end) = u_new(:, end-1); 
    v_new(:, [1, end]) = 0; 

    % No-slip condition at cylinder
    u_new(cylinder_mask) = 0;
    v_new(cylinder_mask) = 0;

    % Pressure correction using central difference
    p_new = p - dt * rho * ((p(:, [2:end, end]) - 2*p + p(:, [1, 1:end-1])) / dx^2 ...
                          + (p([2:end, end], :) - 2*p + p([1, 1:end-1], :)) / dy^2);

    % Update fields
    u = u_new;
    v = v_new;
    p = p_new;

    % Monitor convergence
    if mod(iter, 100) == 0
        fprintf('Iteration %d\n', iter);
    end
end

Drag Coefficient Calculation

We’ll now compute the drag force and determine the drag coefficient:

% Compute drag force components
fx = -trapz(y, p(:, round(Nx/2)));  % Pressure force contribution
tau_wall = mu * (u(:, round(Nx/2)) - 0) / dx; % Shear stress at the wall
fv = trapz(y, tau_wall);  % Viscous force contribution

% Total drag force
Fd = fx + fv;

% Drag coefficient calculation
A = D; % Projected area for 2D cylinder
Cd = Fd / (0.5 * rho * U^2 * A);

fprintf('Drag Coefficient (Cd): %.4f\n', Cd);

Post Processing and Visualization

Finally, we’ll visualize the results with velocity magnitude contours and streamlines.

Velocity Contour Plot

figure;
contourf(X, Y, sqrt(u.^2 + v.^2), 20);
colorbar;
title('Velocity Magnitude Contour');
xlabel('X');
ylabel('Y');
axis equal;

Streamline Plot

figure;
quiver(X, Y, u, v);
title('Velocity Field');
xlabel('X');
ylabel('Y');
axis equal;

Results and Comparison

For low Reynolds numbers, the analytical drag coefficient for a 2D cylinder is:

[ C_d = \frac{8\pi}{\text{Re}} ]

Let’s compare our computed drag coefficient with the analytical value:

Cd_analytical = 8 * pi / Re;
fprintf('Analytical Drag Coefficient: %.4f\n', Cd_analytical);
fprintf('Error: %.2f%%\n', abs(Cd - Cd_analytical) / Cd_analytical * 100);

Further Improvements

To enhance this model, consider:

Refining the grid resolution for better accuracy.
Implementing higher-order numerical schemes.
Extending to 3D simulations.

I hope you found this tutorial helpful. Feel free to try it out and experiment with different parameters!

Drag Coefficent in 2D in Matlab