CFD Online Discussion Forums - 2-D Lid Driven Cavity (Matlab)

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

Code:

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



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

%%

%--------------------------------%

%    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