# 2-D Lid Driven Cavity (Matlab)

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 May 5, 2013, 04:11 2-D Lid Driven Cavity (Matlab) #1 New Member   Join Date: May 2013 Posts: 1 Rep Power: 0 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 04:20. Reason: details

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

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post anna_simons Main CFD Forum 3 August 8, 2012 11:00 new_at_this Main CFD Forum 0 December 22, 2011 17:04 gholamghar Main CFD Forum 0 August 1, 2010 01:55 josephlm Main CFD Forum 3 June 24, 2010 01:58 KK FLUENT 1 December 16, 2009 19:10

All times are GMT -4. The time now is 05:53.

 Contact Us - CFD Online - Top