CFD Online URL
[Sponsors]
Home > Forums > Main CFD Forum

2-D Lid Driven Cavity (Matlab)

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   May 5, 2013, 05:11
Default 2-D Lid Driven Cavity (Matlab)
  #1
New Member
 
Join Date: May 2013
Posts: 1
Rep Power: 0
DaBears13 is on a distinguished road
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

Last edited by DaBears13; May 5, 2013 at 05:20. Reason: details
DaBears13 is offline   Reply With Quote

Reply

Tags
2-d, driven cavity, matlab, navier-stokes solver

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
lid driven cavity matlab code anna_simons Main CFD Forum 3 August 8, 2012 12:00
Lid Driven Cavity velocity profiles new_at_this Main CFD Forum 0 December 22, 2011 17:04
is there any parallel code for the famous Lid Driven Cavity flow? gholamghar Main CFD Forum 0 August 1, 2010 02:55
2D Lid Driven Cavity Flow simulation using MATLAB josephlm Main CFD Forum 3 June 24, 2010 02:58
Lid driven cavity flow KK FLUENT 1 December 16, 2009 19:10


All times are GMT -4. The time now is 02:01.