nx=8;ny=8;dt=0.02;nstep=600;mu=0.1;maxit=100;beta=1.2;h=1/nx; u=zeros(nx,ny+1);v=zeros(nx+1,ny);p=zeros(nx+1,ny+1); ut=zeros(nx,ny+1);vt=zeros(nx+1,ny);c=zeros(nx+1,ny+1)+0.25; uu=zeros(nx,ny);vv=zeros(nx,ny);w=zeros(nx,ny); c(2,3:ny)=1/3;c(nx,3:ny)=1/3;c(3:nx,2)=1/3;c(3:nx,ny)=1/3; c(2,2)=1/2;c(2,ny)=1/2;c(nx,2)=1/2;c(nx,ny)=1/2; un=1;us=0;ve=0;vw=0;time=0.0; for is=1:nstep u(1:nx,1)=2*us-u(1:nx,2);u(1:nx,ny+1)=2*un-u(1:nx,ny); v(1,1:ny)=2*vw-v(2,1:ny);v(nx+1,1:ny)=2*ve-v(nx,1:ny); for i=2:nx-1,for j=2:ny % temporary u-velocity ut(i,j)=u(i,j)+dt*(-(0.25/h)*((u(i+1,j)+u(i,j))^2-(u(i,j)+... u(i-1,j))^2+(u(i,j+1)+u(i,j))*(v(i+1,j)+... v(i,j))-(u(i,j)+u(i,j-1))*(v(i+1,j-1)+v(i,j-1)))+... (mu/h^2)*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)-4*u(i,j))); end,end for i=2:nx,for j=2:ny-1 % temporary v-velocity vt(i,j)=v(i,j)+dt*(-(0.25/h)*((u(i,j+1)+u(i,j))*(v(i+1,j)+... v(i,j))-(u(i-1,j+1)+u(i-1,j))*(v(i,j)+v(i-1,j))+... (v(i,j+1)+v(i,j))^2-(v(i,j)+v(i,j-1))^2)+... (mu/h^2)*(v(i+1,j)+v(i-1,j)+v(i,j+1)+v(i,j-1)-4*v(i,j))); end,end for it=1:maxit % solve for pressure for i=2:nx,for j=2:ny p(i,j)=beta*c(i,j)*(p(i+1,j)+p(i-1,j)+p(i,j+1)+p(i,j-1)-... (h/dt)*(ut(i,j)-ut(i-1,j)+vt(i,j)-vt(i,j-1)))+(1-beta)*p(i,j); end,end end % correct the velocity u(2:nx-1,2:ny)=... ut(2:nx-1,2:ny)-(dt/h)*(p(3:nx,2:ny)-p(2:nx-1,2:ny)); v(2:nx,2:ny-1)=... vt(2:nx,2:ny-1)-(dt/h)*(p(2:nx,3:ny)-p(2:nx,2:ny-1)); time=time+dt % plot the results uu(1:nx,1:ny)=0.5*(u(1:nx,2:ny+1)+u(1:nx,1:ny)); vv(1:nx,1:ny)=0.5*(v(2:nx+1,1:ny)+v(1:nx,1:ny)); w(1:nx,1:ny)=(u(1:nx,2:ny+1)-u(1:nx,1:ny)-... v(2:nx+1,1:ny)+v(1:nx,1:ny))/(2*h); hold off,quiver(flipud(rot90(uu)),flipud(rot90(vv)),'r'); hold on;contour(flipud(rot90(w)),20),axis equal,pause(0.01) end