CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

lid driven cavity

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   July 7, 2007, 08:04
Default 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,

  Reply With Quote

Old   July 8, 2007, 15:55
Default 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.

  Reply With Quote

Old   July 9, 2007, 01:20
Default 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?
  Reply With Quote

Old   July 9, 2007, 12:41
Default 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
  Reply With Quote

Old   July 9, 2007, 12:50
Default 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...
  Reply With Quote

Old   July 11, 2007, 03:05
Default 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?
  Reply With Quote

Old   November 21, 2010, 06:11
Default Solving lid driven cavity problem in MATLAB using Streamfunction-Vorticity method
  #7
New Member
 
Vignesh
Join Date: Nov 2010
Posts: 2
Rep Power: 0
vignesh is on a distinguished road
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
vignesh is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Lid Driven Cavity using Ghost Cell Method and in C++ illuminati5288 Main CFD Forum 0 August 12, 2011 22:05
is there any parallel code for the famous Lid Driven Cavity flow? gholamghar Main CFD Forum 0 August 1, 2010 01:55
2D Lid Driven Cavity Flow simulation using MATLAB josephlm Main CFD Forum 3 June 24, 2010 01:58
Paraview - Lid Driven Cavity kieranhood OpenFOAM Paraview & paraFoam 0 February 13, 2010 17:28
Lid driven cavity flow KK FLUENT 1 December 16, 2009 19:10


All times are GMT -4. The time now is 12:14.