[obdwinston]

I developed a MATLAB code to simulate a simple duct flow, based on the pressure correction technique with staggered grid introduced in J.D. Anderson's textbook on Computational Fluid Dynamics (page 247 - 264).

I managed to get it to work after some trial and error of the space and time steps, relaxation factor, and pressure inflow boundary condition. I am still not exactly sure why it works because I was merely trying my luck with various values. Can someone verify that the code is really doing its job, and not by some fluke that it works.

Thanks in advance!

--- MATLAB CODE ---

clear all, close all, clc

% initialise
m = 25;
n = 25;

p = zeros(m,n);
p_corr = zeros(m,n);
u = zeros(m,n-1);
v = zeros(m-1,n);

% constants
rho = 1.225; % density
nu = 1.5e-5; % kinematic viscosity
dx = 0.0075;
dy = 0.0005;
dt = 0.001;
a_p = .1; % relaxation factor

% boundary conditions
p(1:m,1) = .001; % inflow pressure boundary
p(1:m,2) = .001; % inflow pressure boundary

% SIMPLE algorithm
ctr = 0;
t = 0;
for iter = 1:10000
% x-momentum equation
for i = 2:m-1
for j = 2:n-2
A1 = -((u(i,j+1)*u(i,j+1) - u(i,j-1)*u(i,j-1))/(2*dx) + (u(i-1,j)*.5*(v(i-1,j)+v(i-1,j+1)) - u(i+1,j)*.5*(v(i,j)+v(i,j+1)))/(2*dy));
A2 = nu*((u(i,j+1)-2*u(i,j)+u(i,j-1))/(dx*dx) + (u(i-1,j)-2*u(i,j)+u(i+1,j))/(dy*dy));
A = A1 + A2;
u(i,j) = u(i,j) + A*dt - (dt/dx/rho)*(p(i,j+1)-p(i,j));
end
end

% zeroth-order extrapolation
u(2:m-1,1) = u(2:m-1,2);
u(2:m-1,n-1) = u(2:m-1,n-2);

% y-momentum equation
for i = 2:m-2
for j = 2:n-1
B1 = -((v(i-1,j)*v(i-1,j) - v(i+1,j)*v(i+1,j))/(2*dy) + (v(i,j+1)*.5*(u(i,j)+u(i+1,j)) - v(i,j-1)*.5*(u(i,j-1)+u(i+1,j-1)))/(2*dx));
B2 = nu*((v(i-1,j)-2*v(i,j)+v(i+1,j))/(dy*dy) + (v(i,j+1)-2*v(i,j)+v(i,j-1))/(dx*dx));
B = B1 + B2;
v(i,j) = v(i,j) + B*dt - (dt/dy/rho)*(p(i,j)-p(i+1,j));
end
end

% zeroth-order extrapolation
v(2:m-2,n) = v(2:m-2,n-1);

% pressure correction equation (Poisson equation)
a = 2.*(dt/(dx*dx) + dt/(dy*dy));
b = -dt/(dx*dx);
c = -dt/(dy*dy);

p_corr_diff = ones(m,n);
p_corr_prev = zeros(m,n);
while max(max(abs(p_corr_diff))) > 1e-6
p_corr_prev = p_corr;
for i = 2:m-1
for j = 2:n-1
d = (rho/dx)*(u(i,j) - u(i,j-1)) + (rho/dy)*(v(i-1,j) - v(i,j));
p_corr(i,j) = -(1/a)*(b*p_corr(i,j+1) + b*p_corr(i,j-1) + c*p_corr(i-1,j) + c*p_corr(i+1,j) + d);
end
end
p_corr_diff = p_corr - p_corr_prev;
end

% pressure correction
p = p + a_p*p_corr;
p(1,1:n) = p(2,1:n); % dp/dy = 0 at wall
p(m,1:n) = p(m-1,1:n); % dp/dy = 0 at wall

p(1:m,1) = .1; % inflow pressure boundary
p(1:m,2) = .1; % inflow pressure boundary

ctr = ctr + 1;
t = t + dt;
imagesc(u)
colorbar
str = sprintf('ctr = %d, t = %.3e s', ctr, t);
title(str)
pause(0.001)
end