CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Main CFD Forum (http://www.cfd-online.com/Forums/main/)
-   -   2-D Lid Driven Cavity (Matlab) (http://www.cfd-online.com/Forums/main/117227-2-d-lid-driven-cavity-matlab.html)

 DaBears13 May 5, 2013 04:11

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

 All times are GMT -4. The time now is 10:26.