lid driven cavity
hi all,
I have written a 2D finite volume steady incompressible NavierStokes 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, 
Re: lid driven cavity
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.

Re: lid driven cavity
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? 
Re: lid driven cavity
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 
Re: lid driven cavity
thank you very much harish.. i think it will help me to make my code perfect...

Re: lid driven cavity
hi all,
Now, I have observed that my final solution changes with respect to the underrelaxation factor used.. Also it converges only for very low underrelaxation values (nearly 0.1).. can any one explain me how to handle that factor properly in CFD? 
Solving lid driven cavity problem in MATLAB using StreamfunctionVorticity method
I'm a beginner in CFD and coding; I have to solve the lid driven cavity problem for Re = 150 (nondimensionalised). 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/(M1); DY = 1/(N1); 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 = 1e6; k=0; while(abs(RES)>TOL  k<100), k = k+1; for i=2 : M1, for j=2: N1, R=0.25*(PSIK(i+1,j)+PSIK(i1,j)+PSIK(i,j1)+PSIK(i,j+1)+(ZETAK(i,j)*DX*DX)); PSIK(i,j)=OMEGA*R + (1.0OMEGA)*PSIK(i,j); end end %Calculating velocity field for i=2:M1, for j=2:N1, UK(i,j) = (PSIK(i,j+1)  PSIK(i, j1))/(2*DX); VK(i,j) = (PSIK(i1,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,N1)))); 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((M1),j))+(1.0*(DY)))); end ZETAK1=ZETAK; %Marching forward %Forward euler for i=2:M1, for j=2:N1, ZETAK1(i,j)= H*(((1/RE)*(((ZETAK(i+1,j)2*ZETAK(i,j)+ZETAK(i1,j))/DX^2)+((ZETAK(i,j+1)2*ZETAK(i,j)+ZETAK(i,j1))/DX^2)))(((UK(i,j)*(ZETAK(i+1,j)ZETAK(i1,j)))(VK(i,j)*(ZETAK(i,j+1)ZETAK(i,j1))))/2*DX))+ZETAK(i,j); end end %FTCS method % for i=2:M1, % for j=2:N1, % 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(i1,j)))  ((SIGMAY(i,j)/2)*(ZETAK(i,j+1)ZETAK(i,j1))) + (BETAX*((ZETAK(i+1,j)) (2*ZETAK(i,j)) + (ZETAK(i1,j)))) + (BETAY*((ZETAK(i,j+1)) (2*ZETAK(i,j)) + (ZETAK(i,j1))))); % end % end for i=2 : M1, for j=2: N1, RES=ZETAK1(i,j)*DX*DX + 4*PSIK(i,j)  PSIK(i+1,j)  PSIK(i1,j)  PSIK(i,j+1)  PSIK(i,j1); if (abs(RES)>TOL) ZETAK=ZETAK1; break; end end if (abs(RES)>TOL) break; end end end %Calculating inner point values for i = 2 : M1, for j = 2 : N1, a = 1/4*(2*((PSIK(i+1,j)2*PSIK(i,j)+PSIK(i1,j))*(PSIK(i,j+1)2*PSIK(i,j)+PSIK(i,j1))/DX^2  ((PSIK(i+1,j+1)PSIK(i1,j+1))(PSIK(i+1,j1)+PSIK(i1,j1)))^2/16/DX^2)  P(i+1,j)  P(i1,j)  P(i,j+1)  P(i,j1)); P(i,j) = OMEGA*a + (1.0OMEGA)*P(i,j); if i==2 P(1,j) = 1/3/RE*(ZETAK(1,j+1)ZETAK(1,j1)) + 4/3*P(2,j)  P(3,j)/3; end if i==M1 P(M,j) = 1/3/RE*(ZETAK(M,j+1)ZETAK(M,j1)) + 4/3*P(M1,j)  P(M2,j)/3; end if j==2 P(i,1) = 1/3/RE*(ZETAK(i+1,1)ZETAK(i1,1)) + 4/3*P(i,2)  P(i,3)/3; end if j==M1 P(i,M) = 1/3/RE*(ZETAK(i+1,M)ZETAK(i1,M)) + 4/3*P(i,M1)  P(i,M2)/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:M1,1:N1); 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 
All times are GMT 4. The time now is 23:45. 