CFD Online Discussion Forums - matlab code of marker and cell method for solving 2d incompressible navier stokes eqn

%%MAC algorithm
clc;
close all;
clear all;
%%Initial Constants
l=2; %length of the wind tunnel
w=0.5; %width of the wind tunnel
dx=0.05; %grid size in x direction
dy=0.05; %grid size in y direction
nx=l/dx;
ny=w/dy;
Re=10; %Reynolds number
mu=1/Re; %Coefficient of Viscosity
t=3; %total time
%dt=0.00025; %time step
dt=1;

%Check to ensure that the defined time step meets stability condition?
if dt>0.25*(min(dx,dy))^2/mu || dt>1/Re
dt=min((0.25*(min(dx,dy))^2/mu),(mu));
end
ni=t/dt; %total no of iterations

%Variable initialization-all variables set to zero at initial time
u =zeros(nx+2,ny+2);
utp=zeros(nx+2,ny+2); %Velocity component in x-dir
v =zeros(nx+2,ny+2);
vtp=zeros(nx+2,ny+2); %Velocity component in y-dir
p =zeros(nx+2,ny+2); %Pressure
px=zeros(nx+2,ny+2);
py=zeros(nx+2,ny+2);

%-------------------------------------------------------------------------%
%Time Integration Begins
%-------------------------------------------------------------------------%
for n=1:ni
%Impose Boundary conditions on velocity
%BC's for left wall
u(2,1:ny+2)=1;
u(1,1:ny+2)=2*u(2,1:ny+2)-u(3,1:ny+2);
%BC's for bottom wall
v(1:nx+2,2)=0;
u(1:nx+2,1)=-u(1:nx+2,2);
v(1:nx+2,1)=2*v(1:nx+2,2)-v(1:nx+2,3);
%BC's for Top wall
v(1:nx+2,ny+2)=0;
u(1:nx+2,ny+2)=-u(1:nx+2,ny+1);

%-------------------------------------------------------------------------%
%Advect u and v to get temporary values u_tp and v_tp.
%-------------------------------------------------------------------------%
for i=2:nx+1
for j=2:ny+1

%Pressure at inlet%
p(1,1:ny+2)=101325;
p(2,1:ny+2)=101325;

%First derivatives of velocities%
ux(i,j)= (u(i+1,j)-u(i-1,j))/2*dx;
uy(i,j)= (u(i,j+1)-u(i,j-1))/2*dy;
vx(i,j)= (v(i+1,j)-v(i-1,j))/2*dx;
vy(i,j)= (v(i,j+1)-v(i,j-1))/2*dy;

%Second derivatives of velocities%
u2x(i,j)=(u(i+1,j)-2*u(i,j)+u(i-1,j))/(dx)^2;
u2y(i,j)=(u(i,j+1)-2*u(i,j)+u(i,j-1))/(dy)^2;
v2x(i,j)=(v(i+1,j)-2*v(i,j)+v(i-1,j))/(dx)^2;
v2y(i,j)=(v(i,j+1)-2*v(i,j)+v(i,j-1))/(dy)^2;

%Pressure gradients of cells%
px(i,j)=(p(i,j)-p(i-1,j))/dx;
py(i,j)=(p(i,j)-p(i,j-1))/dy;

%Pressure gradient assumptions at right and bottom boundaries%
px(nx+1,1:ny+2)=px(nx,1:ny+2);
py(1:nx+2,2)=py(1:nx+2,3);

%Temporary velocities in x and y directions%
utp(i,j)=u(i,j)+dt*(-px(i,j)+mu*(u2x(i,j)+u2y(i,j))-(u(i,j)*ux(i,j)+v(i,j)*uy(i,j)));
vtp(i,j)=v(i,j)+dt*(-py(i,j)+mu*(v2x(i,j)+v2y(i,j))-(u(i,j)*vx(i,j)+v(i,j)*vy(i,j)));

end
end

for i=2:nx+1
for j=2:ny+1
D(i,j)=(utp(i,j)-utp(i-1,j))/dx+(vtp(i,j)-vtp(i-1,j))/dy
p1(i,j)=(-D(i,j))/2*dt*((1/dx^2)+(1/dy^2));

while D(i,j)>0.005
utp(i,j)=utp(i,j)+(dt/dx)*p1(i,j);
vtp(i,j)=vtp(i,j)+(dt/dy)*p1(i,j);
utp(i-1,j)=utp(i-1,j)-(dt/dx)*p1(i,j);
vtp(i-1,j)=vtp(i-1,j)-(dt/dy)*p1(i,j);
D(i,j)=(utp(i,j)-utp(i-1,j))/dx+(vtp(i,j)-vtp(i-1,j))/dy;
p1(i,j)=(-D(i,j))/2*dt*((1/dx^2)+(1/dy^2));
end
u(i,j)=utp(i,j);
v(i,j)=vtp(i,j);
end
end

n
end

this is the code i wrote till now.

i am getting very awkward results.
and i think the problem in the formulae i used and i am not able to go forward as i dont know the correct ones. plz help