# 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.  Thread Tools Search this Thread Show Printable Version Email this Page Search this Thread: Advanced Search Display Modes Linear Mode Switch to Hybrid Mode Switch to Threaded Mode 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 OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules Similar Threads Thread Thread Starter Forum Replies Last Post Albert Badal Main CFD Forum 0 February 27, 1999 19:53

All times are GMT -4. The time now is 17:52.