# lid driven cavity

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

 July 7, 2007, 08:04 lid driven cavity #1 sudha Guest   Posts: n/a hi all, I have written a 2D finite volume steady incompressible Navier-Stokes code for lid driven cavity problem. I observed that at Re=100, it gives accurate results for 50*50 grid. But if I increase the number of grid points, it behaves with increased numerical viscosity (max u velocity at the centreline keeps on reducing as I increase the number of grids). I have used central differencing scheme for both convective and diffusive terms. can any one tell me is this the property of the scheme? Also where I can find materials regarding how to implement boundary conditions exactly in staggered grid? thank you,

 July 8, 2007, 15:55 Re: lid driven cavity #2 Harish Guest   Posts: n/a Central difference schemes do not have inbuilt dissipation and hence when you keep refining the grid without adding stability,oscillations can be produced.This was the reason your code worked fine with increased numerical viscosity.You might want to look into the book of patankar ( Numerical heat tranfer and fluid flow) for more information on boundary condition implementation on staggered grid.For stability analysis of numerical schemes any finite difference or CFD books should provide you the information.

 July 9, 2007, 01:20 Re: lid driven cavity #3 sudha Guest   Posts: n/a hi harish, thanks for your reply.. But i cant understand what do you exactly mean by "when you keep refining the grid without adding stability,oscillations can be produced." I read, once we refine the grid, we should get more accurate solution. Can you give me some more explanations about your point?

 July 9, 2007, 12:41 Re: lid driven cavity #4 Harish Guest   Posts: n/a Have a look into any CFD book for some discussion about dissipation and dispersion.It might help you. http://www.cscamm.umd.edu/centpack/publications/#c4

 July 9, 2007, 12:50 Re: lid driven cavity #5 sudha Guest   Posts: n/a thank you very much harish.. i think it will help me to make my code perfect...

 July 11, 2007, 03:05 Re: lid driven cavity #6 sudha Guest   Posts: n/a hi all, Now, I have observed that my final solution changes with respect to the under-relaxation factor used.. Also it converges only for very low under-relaxation values (nearly 0.1).. can any one explain me how to handle that factor properly in CFD?

 November 21, 2010, 06:11 Solving lid driven cavity problem in MATLAB using Streamfunction-Vorticity method #7 New Member   Vignesh Join Date: Nov 2010 Posts: 2 Rep Power: 0 I'm a beginner in CFD and coding; I have to solve the lid driven cavity problem for Re = 150 (non-dimensionalised). I've used the forward euler scheme in the code, though i'm supposed to use FTCS. I tried implementing it but there were some problems and couldn't get the plot properly. Could someone help me to check if the code and plot are correct? And also how to implement FTCS into the code? It would be really helpful if anyone could let me know.. Thank you very much!!! function u = strvor(N) M = N+1; N = N+1; %Setting initial values DX = 1/(M-1); DY = 1/(N-1); RE = 150; OMEGA = 1.9; %Q = 1; %H = min(((RE/2)*(((DX*DX)*(DY*DY))/((DX*DX)+(DY*DY)))), ((2/(RE*(Q*Q))))); H = 0.0005; BETAX = ((1/RE)*(H/(DX*DX))); BETAY = ((1/RE)*(H/(DY*DY))); UK = zeros (M, N); VK = zeros (M, N); P = zeros (M, N); PSIK = zeros (M, N); ZETAK = zeros (M, N); SIGMAX = zeros (M, N); SIGMAY = zeros (M, N); %Setting boundary conditions for i=1:M, UK(i,1) = 0; UK(i,N) = 0; VK(i,1) = 0; VK(i,N) = 0; end for j=1:N, UK(1,j) = 0; UK(M,j) = 1.0; VK(1,j) = 0; VK(M,j) = 0; end %Solving Poisson equation (SOR) RES = 100; TOL = 1e-6; k=0; while(abs(RES)>TOL || k<100), k = k+1; for i=2 : M-1, for j=2: N-1, R=0.25*(PSIK(i+1,j)+PSIK(i-1,j)+PSIK(i,j-1)+PSIK(i,j+1)+(ZETAK(i,j)*DX*DX)); PSIK(i,j)=OMEGA*R + (1.0-OMEGA)*PSIK(i,j); end end %Calculating velocity field for i=2:M-1, for j=2:N-1, UK(i,j) = (PSIK(i,j+1) - PSIK(i, j-1))/(2*DX); VK(i,j) = (PSIK(i-1,j) - PSIK(i+1, j))/(2*DY); end end %Vorticity boundary conditions for i=1:M, ZETAK(i,N) = ((2/(DX*DX))*((PSIK(i,N))-(PSIK(i,N-1)))); ZETAK(i,1) = ((2/(DX*DX))*((PSIK(i,1)-(PSIK(i,2))))); end for j=1:N, ZETAK(1,j) = ((2/(DX*DX))*((PSIK(1,j)-(PSIK(2,j))))); ZETAK(M,j) = ((2/(DY*DY))*((PSIK(M,j))-(PSIK((M-1),j))+(1.0*(-DY)))); end ZETAK1=ZETAK; %Marching forward %Forward euler for i=2:M-1, for j=2:N-1, ZETAK1(i,j)= H*(((1/RE)*(((ZETAK(i+1,j)-2*ZETAK(i,j)+ZETAK(i-1,j))/DX^2)+((ZETAK(i,j+1)-2*ZETAK(i,j)+ZETAK(i,j-1))/DX^2)))-(((UK(i,j)*(ZETAK(i+1,j)-ZETAK(i-1,j)))-(VK(i,j)*(ZETAK(i,j+1)-ZETAK(i,j-1))))/2*DX))+ZETAK(i,j); end end %FTCS method % for i=2:M-1, % for j=2:N-1, % SIGMAX(i,j) = ((UK(i,j) * H)/DX); % SIGMAY(i,j) = ((VK(i,j) * H)/DY); % ZETAK1(i,j) = ((ZETAK(i,j)) - ((SIGMAX(i,j)/2)*(ZETAK(i+1,j)-ZETAK(i-1,j))) - ((SIGMAY(i,j)/2)*(ZETAK(i,j+1)-ZETAK(i,j-1))) + (BETAX*((ZETAK(i+1,j))- (2*ZETAK(i,j)) + (ZETAK(i-1,j)))) + (BETAY*((ZETAK(i,j+1))- (2*ZETAK(i,j)) + (ZETAK(i,j-1))))); % end % end for i=2 : M-1, for j=2: N-1, RES=-ZETAK1(i,j)*DX*DX + 4*PSIK(i,j) - PSIK(i+1,j) - PSIK(i-1,j) - PSIK(i,j+1) - PSIK(i,j-1); if (abs(RES)>TOL) ZETAK=ZETAK1; break; end end if (abs(RES)>TOL) break; end end end %Calculating inner point values for i = 2 : M-1, for j = 2 : N-1, a = -1/4*(2*((PSIK(i+1,j)-2*PSIK(i,j)+PSIK(i-1,j))*(PSIK(i,j+1)-2*PSIK(i,j)+PSIK(i,j-1))/DX^2 - ((PSIK(i+1,j+1)-PSIK(i-1,j+1))-(PSIK(i+1,j-1)+PSIK(i-1,j-1)))^2/16/DX^2) - P(i+1,j) - P(i-1,j) - P(i,j+1) - P(i,j-1)); P(i,j) = OMEGA*a + (1.0-OMEGA)*P(i,j); if i==2 P(1,j) = 1/3/RE*(ZETAK(1,j+1)-ZETAK(1,j-1)) + 4/3*P(2,j) - P(3,j)/3; end if i==M-1 P(M,j) = 1/3/RE*(ZETAK(M,j+1)-ZETAK(M,j-1)) + 4/3*P(M-1,j) - P(M-2,j)/3; end if j==2 P(i,1) = -1/3/RE*(ZETAK(i+1,1)-ZETAK(i-1,1)) + 4/3*P(i,2) - P(i,3)/3; end if j==M-1 P(i,M) = -1/3/RE*(ZETAK(i+1,M)-ZETAK(i-1,M)) + 4/3*P(i,M-1) - P(i,M-2)/3; end end end %Plotting figure(1)% Streamline plot with number Z = PSIK(1:M,1:N); X = linspace(0,1,size(Z,2)); Y = linspace(0,1,size(Z,1)); [c,h] = contour(X,Y,Z); figure(2)% Pressure plot Z = P(1:M-1,1:N-1); X = linspace(0,1,size(Z,2)); Y = linspace(0,1,size(Z,1)); [c,h] = contourf(X,Y,Z); axis equal axis([0 1 0 1]) drawnow u=VK; end

 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 illuminati5288 Main CFD Forum 0 August 12, 2011 22:05 gholamghar Main CFD Forum 0 August 1, 2010 01:55 josephlm Main CFD Forum 3 June 24, 2010 01:58 kieranhood OpenFOAM Paraview & paraFoam 0 February 13, 2010 17:28 KK FLUENT 1 December 16, 2009 19:10

All times are GMT -4. The time now is 21:28.

 Contact Us - CFD Online - Top