CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

matlab code of marker and cell method for solving 2d incompressible navier stokes eqn

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   November 10, 2011, 09:55
Default matlab code of marker and cell method for solving 2d incompressible navier stokes eqn
  #1
New Member
 
sandeep gangarapu
Join Date: Nov 2011
Posts: 1
Rep Power: 0
sgangarapu is on a distinguished road
%%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
sgangarapu is offline   Reply With Quote

Old   August 28, 2013, 10:33
Smile Inquiry
  #2
New Member
 
Join Date: Aug 2013
Posts: 1
Rep Power: 0
choclatebar is on a distinguished road
Hi,
I am currently doing a similar project. Have you solved the problem? What's the meaning of p1? I am looking forward to your reply. Thank you.
choclatebar is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Marker And Cell method (MAC) Albert Badal Main CFD Forum 0 February 27, 1999 20:53


All times are GMT -4. The time now is 02:15.