# 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

 November 10, 2011, 08:55 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 %%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

 August 28, 2013, 10:33 Inquiry #2 New Member   Join Date: Aug 2013 Posts: 1 Rep Power: 0 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.