CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Main CFD Forum (https://www.cfd-online.com/Forums/main/)
-   -   lid driven cavity (https://www.cfd-online.com/Forums/main/13779-lid-driven-cavity.html)

 sudha July 7, 2007 08:04

lid driven cavity

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,

 Harish July 8, 2007 15:55

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.

 sudha July 9, 2007 01:20

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?

 Harish July 9, 2007 12:41

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

 sudha July 9, 2007 12:50

Re: lid driven cavity

thank you very much harish.. i think it will help me to make my code perfect...

 sudha July 11, 2007 03:05

Re: lid driven cavity

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?

 vignesh November 21, 2010 06:11

Solving lid driven cavity problem in MATLAB using Streamfunction-Vorticity method

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

 All times are GMT -4. The time now is 22:45.