[DaBears13]

Hi,
I have to write a Navier-Stokes solver for a 2-D Lid Driven Cavity. My teacher gave me a portion of his code (the Poisson pressure solver and some I.C.'s and B.C.'s) and I have to fill in the blanks. I've been using http://math.mit.edu/cse/codes/mit18086_navierstokes.pdf, http://sitemaker.umich.edu/anand/files/project_1.pdf, and http://www.mathworks.com/matlabcentr...en-cavity-flow to help me out but I can't get the visualization right. I want to have a quiver plot of the velocity vectors, streamline plot, and contour plot of the pressure. Any help is appreciated! Thanks in advance!

Below is my code

Edit* I am using a FTCS MAC method on a staggered grid

%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%
%%
%--------------------------------%
%    Solve the NS Equations      %
%       using MAC Method         %
%--------------------------------%

%%%%%%%%%% Setup %%%%%%%%%%
close all;
clear all;

% input
N          = 100;    % number of grid point
Final_Time = 20.0;   % final time
dtout      = 1.0;
Re         = 400;    % Reynolds number
cfl        = 1.0;    % cfl criterion

% allocate variables

p = ones(N,N);
u = zeros(N+1,N);
v = zeros(N,N+1);
F = zeros(N+1,N);
G = zeros(N,N+1);
Q = zeros(N-2,N-2);

% grid
L     = 1;
dx    = L/(N-1);
dy    = L/(N-1);
x     = 0:dx:L;
y     = 0:dy:L;

% time step
dt        = 0.01;                % I chose
Nsteps    = floor(Final_Time/dt);
nout      = floor(dtout/dt);    % number of time steps for output

% Poisson matrix
rf    = 1.6;
eps   = 0.001;
itmax = 5000;

% velocity profiles
UN = 1;     VN = 0;
US = 0;     VS = 0;
UE = 0;     VE = 0;
UW = 0;     VW = 0;

%%
%%%%%%%%%% main time loop %%%%%%%%%%
for tstep = 1:Nsteps
 
  % boundary condtions %
  
  % u-vel
    u(1,:)      = 2*UW - u(2,:);    % Left wall
    u(end,:)    = 2*UE - u(end-1,:);% Right wall
    u(:,1)      = US;               % Bottom wall
    u(:,end)    = UN;               % Top wall
    
  % v-vel
    v(1,:)      = VW;               % Left wall
    v(end,:)    = VE;               % Right wall
    v(:,1)      = 2*VS - v(:,2);    % Bottom wall
    v(:,end)    = 2*VN - v(:,end-1);% Top wall
    
    
%%%%%%%%%% compute f,g,q %%%%%%%%%%

%%%%% f

    Conv_u = u; Diff_u = u;
    
    i = 2:N-1;
    j = 2:N-1;

    Diff_u(i+1,j)   = ((u(i+2,j) - 2.*u(i+1,j) + u(i,j))/dx^2) ...
                        + ((u(i+1,j+1) - 2.*u(i+1,j) + u(i+1,j-1))/dy^2);
    
    Conv_u(i+1,j)   = ((((u(i+2,j)+u(i+1,j))/2).^2-((u(i+1,j)+u(i,j))/2).^2)/dx) ...
                        + ((((u(i+1,j)+u(i+1,j+1))/2).*((v(i,j+1)+v(i+1,j+1))/2) ...
                        - ((u(i+1,j)+u(i+1,j-1))/2).*((v(i,j)+v(i+1,j))/2))/dy);
    
    F(i+1,j)        = u(i+1,j) + dt*(Diff_u(i+1,j)/Re - Conv_u(i+1,j));
    
%%%%% g

    Conv_v = v; Diff_v = v;    

    Diff_v(i,j+1)   = ((v(i+1,j+1) - 2*v(i,j+1) + v(i-1,j+1))/dx^2) ...
                        + ((v(i,j+2) - 2*v(i,j+1) + v(i,j))/dy^2);
    
    Conv_v(i,j+1)   = ((((v(i,j+2)+v(i,j+1))/2).^2-((v(i,j+1)+v(i,j))/2).^2)/dy) ...
                        + ((((v(i,j+1)+v(i+1,j+1))/2).*((u(i+1,j)+u(i+1,j+1))/2) ...
                        - ((v(i,j+1)+v(i-1,j+1))/2).*((u(i,j)+u(i,j-1))/2))/dx);
                    
    G(i,j+1)        = v(i,j+1) + dt*(Diff_v(i,j+1)/Re - Conv_v(i,j+1));
    
    
            
%%%%% q

    Q(i-1,j-1)       = (1/dt) * ((F(i+1,j) - F(i,j))/dx + (G(i,j+1) ...
                        - G(i,j))/dy);

    % Poisson solve
    pold   = p;
    change = 2*eps;
    it     = 1;
    while (change > eps)
      
      pold = p;
    % boundary condition
      p(2:N-1, 1) = p(2:N-1, 2) - 2.0*v(2:N-1, 3)/(Re*dx); % bottom
      p(2:N-1, N) = p(2:N-1,N-1) + 2*v(2:N-1,N-2)/(Re*dx); % top
      p(1, 2:N-1) = p(2, 2:N-1) - 2.0*u(3, 2:N-1)/(Re*dy); % left
      p(N, 2:N-1) = p(N-1,2:N-1) + 2*u(N-2,2:N-1)/(Re*dy); % right
      
      % SOR
        for i=2:N-1
            for j=2:N-1
                p(i,j) = 0.25*(p(i-1,j)+p(i,j-1)+p(i+1,j)+p(i,j+1) - ...
                            Q(i-1,j-1)*dx^2);
                p(i,j) = pold(i,j) + rf*(p(i,j)-pold(i,j));
            end
        end
    
      
        pmax = max(abs(pold(:)));
        if (pmax == 0) 
            pmax = 1.0;
        end
        change =  max(abs( (pold(:)- p(:))/pmax ));
        it = it + 1;
        if (it > itmax) 
            break;
        end    
    end

%%%%%%%%%% output %%%%%%%%%%

    if (mod(tstep,nout) == 0)
        time = tstep*dt;
        % plot output
        uplot(1:N,1:N)  = (1/2) * (u(1:N,1:N) + u(2:N+1,1:N));
        vplot(1:N,1:N)  = (1/2) * (v(1:N,1:N) + v(1:N,2:N+1));
        Len             = sqrt(uplot.^2 + vplot.^2 + eps);
        uplot           = uplot./Len;
        vplot           = vplot./Len;
        q               = ones(N,N);
        q(2:N-1,2:N-1)  = reshape(Q,N-2,N-2);
        % Quiver Plot 
        figure(1)
%         contour(x,y,q',200,'k-'), hold on
%         contour(x,y,Q',50,'k-'), hold on
%         quiver(x,y,uplot',vplot',10,'k-')
%         quiver(x,y,uplot',vplot')
%         streamline(x,y,uplot',vplot',
        sx = 0:.05:2;
        sy = 0:.05:2;
        fn = stream2(x,y,uplot',vplot',sx,sy);
        clf, streamline(fn);
        axis equal
        axis([0 L 0 L])
%         contourf(x,y,p',20,'w-');
        hold off
        colormap(jet)
        colorbar
%         axis equal
        axis([0 L 0 L])
%         p = sort(p); caxis(p([8 end-7]))
        title({['2D Cavity Flow with Re = ',num2str(Re)];
            ['time(\itt) = ',num2str(time)]})
        xlabel('Spatial Coordinate (x) \rightarrow')
        ylabel('Spatial Coordinate (y) \rightarrow')
        drawnow;
        
    end
    
%%%%%%%%%% update velocity %%%%%%%%%%

%%%%% u velocity
    i = 2:N-2;
    j = 2:N-1;
    
    u(i+1,j)    = F(i+1,j) - (dt/dx)*(p(i+1,j) - p(i,j));

%%%%% v velocity
    
    v(i,j+1)    = G(i,j+1) - (dt/dy)*(p(i,j+1) - p(i,j));
    
    
end