1.
Problem Setup and Grid
This MATLAB code solves a 2D incompressible flow problem using the vorticity–stream function
formulation of Navier–Stokes equations. The domain is a square (1 × 1) and is divided into a
structured grid using the Finite Difference Method (FDM). Each point in the grid is used to
approximate derivatives numerically instead of using continuous equations.
2. Flow Variables
Three main variables are used:
Stream function (ψ): defines the flow pattern (streamlines)
Vorticity (ω): represents local rotation of fluid
Velocity (u, v): describes fluid motion in x and y directions
All variables are initially set to zero, and the solution is improved iteratively.
3. Poisson Equation (Stream Function)
The stream function is obtained by solving the Poisson equation, which links ψ and ω. This is solved
using the Gauss–Seidel method, where each grid point is updated using neighbouring values and
current vorticity.
In simple terms:
This step reconstructs the flow pattern from the rotation of fluid.
4. Velocity Calculation
Once ψ is known, velocity is calculated using finite differences:
u is obtained from change of ψ in y-direction
v is obtained from change of ψ in x-direction
This converts the flow representation into actual fluid motion.
5. Boundary Conditions
A lid-driven cavity is used:
Top wall moves with constant velocity
Other walls are stationary (no-slip condition)
This creates circulation and vortex formation inside the domain.
6. Vorticity Transport Equation
Vorticity is updated using the vorticity transport equation, which includes:
Convection (transport by flow)
Diffusion (viscous spreading)
Finite difference approximations are used for spatial derivatives, and explicit time stepping updates the
solution.
7. Residual Verification
After each iteration, the code compares the updated stream function (ψ) with its previous value (ψ old).
The residual is calculated as the maximum absolute difference between these two values over the whole
grid:
Residual = max |ψ(new) − ψ(old)|
If this residual becomes smaller than 10⁻⁶, the iteration stops; otherwise, the loop continues. This check
is applied at the end of each iteration to ensure that the solution of the Poisson equation is converging
properly.
8. Results Interpretation
The code generates:
Stream function contours (flow pattern)
Streamlines (fluid motion paths)
Vorticity contours (rotation strength)
Velocity field (flow direction)
Residual plot (convergence behaviour)
These results help visualize how fluid moves and stabilizes inside the domain.
MATLAB CODE:
clc
clear
close all
%% Domain
nx=81;
ny=81;
L=1;
dx=L/(nx-1);
dy=dx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
[X,Y]=meshgrid(x,y);
%% Parameters
Re=100;
U=1;
dt=0.001;
nt=5000;
tol=1e-6;
%% Variables
psi=zeros(ny,nx);
omega=zeros(ny,nx);
u=zeros(ny,nx);
v=zeros(ny,nx);
residual=zeros(nt,1);
%% Main Loop
for n=1:nt
psi_old_global=psi;
%% Poisson solver (more iterations)
for k=1:200
psi_old=psi;
for i=2:ny-1
for j=2:nx-1
psi(i,j)=0.25*( ...
psi(i+1,j)+psi(i-1,j)+ ...
psi(i,j+1)+psi(i,j-1)+ ...
dx^2*omega(i,j));
end
end
end
%% Velocity from streamfunction
[dpsidy,dpsidx]=gradient(psi,dy,dx);
u=dpsidy;
v=-dpsidx;
%% Lid velocity BC
u(end,:)=U;
u(1,:)=0;
u(:,1)=0;
u(:,end)=0;
v(end,:)=0;
v(1,:)=0;
v(:,1)=0;
v(:,end)=0;
%% Vorticity transport
omega_old=omega;
[dwdy,dwdx]=gradient(omega_old,dy,dx);
convection = u.*dwdx + v.*dwdy;
diffusion = del2(omega_old,dx,dy);
omega = omega_old + dt*( ...
-convection + (1/Re)*diffusion );
%% Vorticity BC
omega(:,1)= -2*psi(:,2)/dx^2;
omega(:,end)= -2*psi(:,end-1)/dx^2;
omega(1,:)= -2*psi(2,:)/dy^2;
omega(end,:)= -2*psi(end-1,:)/dy^2 ...
-2*U/dy;
%% Residual
residual(n)=max(max(abs(psi-psi_old_global)));
if residual(n)<tol && n>500
residual=residual(1:n);
break
end
end
%% -------- PLOTS --------
figure
contourf(X,Y,psi,40)
colorbar
xlabel('x (m)')
ylabel('y (m)')
title('Stream Function Contours')
axis equal tight
figure
contourf(X,Y,omega,40)
colorbar
xlabel('x (m)')
ylabel('y (m)')
title('Vorticity Contours')
axis equal tight
figure
streamslice(X,Y,u,v)
xlabel('x (m)')
ylabel('y (m)')
title('Flow Streamlines')
axis equal tight
figure
quiver(X,Y,u,v)
xlabel('x (m)')
ylabel('y (m)')
title('Velocity Vector Field')
axis equal tight
figure
semilogy(residual,'LineWidth',2)
grid on
xlabel('Iteration')
ylabel('Residual')
title('Residual Convergence')