MATLAB code: Finite difference projection method for 2D incompressible Navier-Stokes equations
Date:
%The space is discretized with finite difference method. The first-order pressure projection algorithm is used to update velocities and pressure. The intermediate velocities are updated with an explicit Euler method. The pressure Poisson equation is solved with a Gauss-Seidel-type iteration.
clear all; clc; close all;
nx = 64; ny = nx;
xleft = 0; xright = 1; yleft = 0; yright = ny/nx*xright;
h = (xright - xleft)/nx; h2 = h^2;
bdvel = 1; % boundary velocity at top
x = linspace(xleft+0.5h, xright-0.5h, nx);
y = linspace(yleft+0.5h, yright-0.5h, ny);
max_it = 10000; ns = max_it/10;
Re = 10000; dt = 0.01h2Re; tol = 0.00001;
% initialization
p = zeros(nx,ny); u = zeros(nx+1,ny+2); v = zeros(nx+2,ny+1);
nu=u; nv=v; tu = zeros(nx+1,ny); tv = zeros(nx,ny+1);
adv_u = tu; adv_v = tv; uu = zeros(nx,ny); vv = zeros(nx,ny);
f = zeros(nx,ny);
for it = 1:max_it
% boundary condition of u
for i = 1:nx+1
u(i,1) = -u(i,2); u(i,ny+2) = 2*bdvel - u(i,ny+1);
end
% boundary condition of v
for j = 1:ny+1
v(1,j) = -v(2,j); v(nx+2,j) = -v(nx+1,j);
end
% advection term using the upwind scheme
for i = 2:nx
for j = 2:ny+1
if u(i,j) > 0
adv_u(i,j-1) = u(i,j)*(u(i,j) - u(i-1,j))/h;
else
adv_u(i,j-1) = u(i,j)*(u(i+1,j) - u(i,j))/h;
end
vc = 0.25*(v(i,j-1) + v(i+1,j-1) + v(i,j) + v(i+1,j));
if vc > 0
adv_u(i,j-1) = adv_u(i,j-1) + vc*(u(i,j) - u(i,j-1))/h;
else
adv_u(i,j-1) = adv_u(i,j-1) + vc*(u(i,j+1) - u(i,j))/h;
end
end
end
for i = 2:nx+1
for j = 2:ny
if v(i,j) > 0
adv_v(i-1,j) = v(i,j)*(v(i,j) - v(i,j-1))/h;
else
adv_v(i-1,j) = v(i,j)*(v(i,j+1) - v(i,j))/h;
end
uc = 0.25*(u(i-1,j) + u(i-1,j+1) + u(i,j) + u(i,j+1));
if uc > 0
adv_v(i-1,j) = adv_v(i-1,j) + uc*(v(i,j) - v(i,j-1))/h;
else
adv_v(i-1,j) = adv_v(i-1,j) + uc*(v(i+1,j) - v(i,j))/h;
end
end
end
% temporal velocity u
for i = 2:nx
for j = 2:ny+1
tu(i,j-1)=u(i,j)+dt*(-adv_u(i,j-1)+(u(i+1,j) ...
+ u(i-1,j) - 4*u(i,j) + u(i,j+1) + u(i,j-1))/(Re*h2));
end
end
% temporal velocity v
for i = 2:nx+1
for j = 2:ny
tv(i-1,j)=v(i,j)+dt*(-adv_v(i-1,j)+(v(i+1,j)+v(i-1,j) ...
-4*v(i,j)+v(i,j+1)+v(i,j-1))/(Re*h2));
end
end
% Poisson equation / source term
for i = 1:nx
for j = 1:ny
f(i,j)=(tu(i+1,j)-tu(i,j)+tv(i,j+1)-tv(i,j))/(h*dt);
end
end
err = 1.0; pold = p; count = 0;
% relaxation using the Gauss-Seidel iterative method
while ( err > tol )
count = count + 1;
for i = 1:nx
for j = 1:ny
sor = f(i,j);
if i == 1
sor = sor - p(i+1,j)/h2; cof = -1;
elseif i == nx
sor = sor - p(i-1,j)/h2; cof = -1;
else
sor = sor - (p(i+1,j) + p(i-1,j))/h2; cof = -2;
end
if j == 1
sor = sor - p(i,j+1)/h2; cof = cof - 1;
elseif j == ny
sor = sor - p(i,j-1)/h2; cof = cof - 1;
else
sor = sor - (p(i,j+1) + p(i,j-1))/h2; cof = cof - 2;
end
p(i,j) = sor/cof*h2;
end
end
% update pressure
p = p - sum(sum(p))/(nx*ny);
err = norm(pold-p,2);
pold = p;
end
count
% updating u
for i = 2:nx
for j = 2:ny+1
nu(i,j)=tu(i,j-1)-dt*(p(i,j-1)-p(i-1,j-1))/h;
end
end
% updating v
for i = 2:nx+1
for j = 2:ny
nv(i,j)=tv(i-1,j)-dt*(p(i-1,j)-p(i-1,j-1))/h;
end
end
u = nu; v = nv;
it
% post-processing
if mod(it,ns) == 0
figure(it); clf; hold on
for i = 2:nx+1
for j = 2:ny+1
uu(i-1,j-1) = 0.5*(u(i-1,j) + u(i,j));
vv(i-1,j-1) = 0.5*(v(i,j-1) + v(i,j));
end
end
sh = 2; s = 0.2;
q = quiver(x(1:sh:end),y(1:sh:end), ...
s*uu(1:sh:end,1:sh:end)',s*vv(1:sh:end,1:sh:end)',0);
st=streamline(x,y,uu',vv',rand(10,1),rand(10,1),[0.1 500]);
set(st,'color','red','linewidth',1);
axis image; box on; axis([0 1 0 1]);
end end