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

lid-driven cavity using Matlab

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

Like Tree2Likes
  • 2 Post By idris

Reply
 
LinkBack Thread Tools Display Modes
Old   November 27, 2006, 10:49
Default lid-driven cavity using Matlab
  #1
M.A
Guest
 
Posts: n/a
Does anyone happen to have a Matlab code for solving the 2-D lid-driven cavity problem in primitive variables (i.e. pressure and velocity)? I will be so grateful, if you would be able to share it.

Thanks a lot.

  Reply With Quote

Old   November 27, 2006, 11:38
Default Re: lid-driven cavity using Matlab
  #2
Andrew Hayes
Guest
 
Posts: n/a
does it have to be in primitive variables? It is not difficult at all using streamfuntion-vorticity, or streamfunction-velocity equations.
  Reply With Quote

Old   November 27, 2006, 12:06
Default Re: lid-driven cavity using Matlab
  #3
M.A
Guest
 
Posts: n/a
Thanks Andrew. In fact, I have already solved the problem using vorticity-streamfuntion formulation (which is quite easy), but I'm having some troubles implementing the SIMPLE algorithm (on staggered grid) for this problem in Matlab.

  Reply With Quote

Old   November 28, 2006, 08:46
Default Re: lid-driven cavity using Matlab
  #4
Andrew Hayes
Guest
 
Posts: n/a
I have you tried using the SIMPLER method? Just out of curiousity, how are you defining the pressure in the cavity? I was reading that the SIMPLE method cannot be used in all situations. Can it be used here? You might need to incorporate the under-relaxation constant.
  Reply With Quote

Old   August 24, 2009, 02:51
Default
  #5
New Member
 
Join Date: Aug 2009
Posts: 5
Rep Power: 8
idris is on a distinguished road
Send a message via Yahoo to idris
Quote:
Originally Posted by Andrew Hayes
;46703
does it have to be in primitive variables? It is not difficult at all using streamfuntion-vorticity, or streamfunction-velocity equations.
hi there, i know those equation, i applied the central difference method and put it in matlab code. The stream function i design it to be dimensionless and the same for the vorticity.

the problem is, after the coding, i can't get the nice solution.

i'm sorry to bother you all, but can you send me the code that's working because my code isn't. I just want to compare or to determine where is my code (matlab or fortran) went wrong. ok.

my email: sesenangje@yahoo.com

Thanx very much
idris is offline   Reply With Quote

Old   August 24, 2009, 09:15
Default
  #6
Member
 
Jonas T. Holdeman, Jr.
Join Date: Mar 2009
Location: Knoxville, Tennessee
Posts: 94
Rep Power: 9
Jonas Holdeman is on a distinguished road
I have a Hermite FE code used in a paper to appear in IJNMF. It computes the velocity and the pressure as a function of the velocity.
Jonas Holdeman is offline   Reply With Quote

Old   August 26, 2009, 00:48
Default
  #7
New Member
 
Join Date: Aug 2009
Posts: 5
Rep Power: 8
idris is on a distinguished road
Send a message via Yahoo to idris
Quote:
Originally Posted by Jonas Holdeman View Post
I have a Hermite FE code used in a paper to appear in IJNMF. It computes the velocity and the pressure as a function of the velocity.
thanx very much, but mine in FD not FE. but if you want to send me the codes, it will be really nice of you. With the codes, i can validate my streamline of a lid driven cavity by comparing your solution with mine. In addition, your solution was accepted by the international journal which i can rely on. thanks for your help. Would you mind send the FE code to me? my email is as shown above

Thank you very very much
idris is offline   Reply With Quote

Old   August 26, 2009, 13:02
Default
  #8
Member
 
Jonas T. Holdeman, Jr.
Join Date: Mar 2009
Location: Knoxville, Tennessee
Posts: 94
Rep Power: 9
Jonas Holdeman is on a distinguished road
Idris,

I will send the Matlab LDC code for one of the methods to you as attachments. A preprint of the paper describing the method can be found at http://webpages.charter.net/jtholdeman/RevIJNMFart1.pdf .
Jonas Holdeman is offline   Reply With Quote

Old   September 2, 2009, 04:30
Default
  #9
New Member
 
Join Date: Aug 2009
Posts: 5
Rep Power: 8
idris is on a distinguished road
Send a message via Yahoo to idris
thanx guys,
i finally created the code myself. it's the matter of boundary condition.
if the boundary condition is wrong, the solution also wrong. Actually my task is to solve vorticity-streamfunction equation of a lid driven cavity at steady state(using central finite difference). and not to forget,the relaxation scheme is needed when the reynolds number is very high (Re >1000)

Anyway, many thanks to everybody, especially Jonas Holdeman.

bye.
idris is offline   Reply With Quote

Old   October 8, 2009, 06:16
Default
  #10
New Member
 
Join Date: Oct 2009
Posts: 7
Rep Power: 7
bhagat is on a distinguished road
can u please share the code....

thank you.
bhagat is offline   Reply With Quote

Old   October 8, 2009, 23:28
Default
  #11
New Member
 
Join Date: Aug 2009
Posts: 5
Rep Power: 8
idris is on a distinguished road
Send a message via Yahoo to idris
ok, here it is

clear all;clf;clc;
re=100; % for higher reynolds number(>500-1000), use under relaxation
rex=1; % underrelaxation, 0<rex<1
rexb=rex;

nx=51; dx=1/(nx-1); %100x100 grid (101*101 nodes)
ny=51; dy=1/(ny-1);

fac=2*(1/dx^2+1/dy^2);
Y=1;
uo=1;
nu=Y/(uo*re);


% initial condition

p(1:nx,1:ny)=0;
o(1:nx,1:ny)=0;
psi(1:nx,1:ny)=0;
omega(1:nx,1:ny)=0;
uact(1:nx,1:ny)=0;
vact(1:nx,1:ny)=0;
u(1:nx,1:ny)=0;
v(1:nx,1:ny)=0;
u(1:nx,ny)=1;

for iteration =1:2000%iteration

disp(['iteration= ',int2str(iteration)])
% stream function(internal nodes)

for i=2nx-1)
for j=2ny-1)
pold(i,j)=p(i,j);
f1=1/fac*(o(i,j)+(p(i+1,j)+p(i-1,j))/dx^2+(p(i,j+1) ...
+p(i,j-1))/dy^2);
p(i,j)=p(i,j)+rex*(f1-p(i,j));
end
end

% Vorticity (Boundary nodes)
for j=1:ny
fab=-2*p(2,j)/dx^2; %left side(AB)
o(1,j)=o(1,j)+rexb*(fab-o(1,j));
fcd=-2*p(nx-1,j)/dx^2; %right side(CD)
o(nx,j)=o(nx,j)+rexb*(fcd-o(nx,j));
end

for i=1:nx
fad=-2*p(i,2)/dy^2; %bottom side(AD)
o(i,1)=o(i,1)+rexb*(fad-o(i,1));
fbc=-(2*p(i,ny-1)+2*uo*dy)/dy^2; %top side (BC)
o(i,ny)=o(i,ny)+rexb*(fbc-o(i,ny));
end

%vorticity internal nodes

for i=2nx-1)
for j=2ny-1)
oold(i,j)=o(i,j);
f2=1/fac*((o(i+1,j)+o(i-1,j))/dx^2+(o(i,j+1)+o(i,j-1))/dy^2 ...
-re*(p(i,j+1)-p(i,j-1))*(o(i+1,j)-o(i-1,j))/(4*dx*dy)...
+re*(p(i+1,j)-p(i-1,j))*(o(i,j+1)-o(i,j-1))/(4*dx*dy));

o(i,j)=o(i,j)+rex*(f2-o(i,j));


end
end

% Calculate the velocity (u and v)

for i=2nx-1)
for j=2ny-1)
u(i,j)=(p(i,j+1)-p(i,j-1))/(2*dy);
v(i,j)=-(p(i+1,j)-p(i-1,j))/(2*dx);
end
end

% rearrange in cartesian coordinate
%(which originally in matrix coordinate)
for i=1:nx
for j=1:ny
k=i;
l=j;
psi(l,k)=p(i,j);
omega(l,k)=o(i,j);
uact(l,k)=u(i,j);
vact(l,k)=v(i,j);
end
end
figure(1)% Streamline plot with number
Z = psi(1:ny,1:nx);
X = linspace(0,1,size(Z,2));
Y = linspace(0,1,size(Z,1));
[c,h] = contour(X,Y,Z);

axis equal
axis([0 1 0 1])
drawnow
end

hope you can further developed these code by yourself. It's up to you what you are going to do about it :-)
metmet and ananthulasriram like this.
idris is offline   Reply With Quote

Old   October 9, 2009, 04:45
Default
  #12
New Member
 
Join Date: Oct 2009
Posts: 7
Rep Power: 7
bhagat is on a distinguished road
thank you very much
bhagat is offline   Reply With Quote

Old   November 13, 2010, 05:55
Default
  #13
Senior Member
 
KHB
Join Date: Aug 2010
Location: Singapore
Posts: 107
Rep Power: 7
lava12005 is on a distinguished road
Hi guys, I am a newbie in CFD thing. Currently I am trying to solve the same problem also.

I solved the problem using the stream function formulation, and now I'm trying to solve in in terms of primitive variable (u,v, and p) using projection method.

First I tried using the most simple explicit projection method and it turns out the solution diverge, and I also used implicit projection method where the value of velocity field in each time step is iterated first and it still diverge. I don't know whats wrong with the code.. Have any of you tried to solve the problem using Projection method?
lava12005 is offline   Reply With Quote

Old   November 15, 2010, 15:29
Default
  #14
New Member
 
Marshall Ericson
Join Date: Nov 2010
Posts: 1
Rep Power: 0
Marshall is on a distinguished road
thank you very much for posting these codes. but these codes seem to have some problems. it doesnt give correct plot nor correct data. if you have any other simple matlab code for lid driven cavity problem that worked , i would be very grateful.
Marshall is offline   Reply With Quote

Old   November 16, 2010, 13:05
Default
  #15
New Member
 
Tanay Deshpande
Join Date: Aug 2010
Posts: 20
Rep Power: 7
Tanay is on a distinguished road
Hello All,
I'm trying to implement SIMPLE in Finite Volume method in driven-cavity flow using C++ in basic variables. Since the problem is steady state, I start out with u=0, v=0, p=0 for all interior nodes. This is causing the code to diverge away to infinity. For high viscosities, the divergence is less pronounced but u doesn't ever become negative in the results as it should. What may be the problem? Is it because I'm using unrelaxed pressure and velocity corrector equations?
Thanks in advance.
Tanay is offline   Reply With Quote

Old   November 20, 2010, 17:40
Default Lid driven cavity problem with Stream-function vorticity formulation for Re = 150
  #16
New Member
 
Vignesh
Join Date: Nov 2010
Posts: 2
Rep Power: 0
vignesh is on a distinguished road
Hello! 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.. 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;



UK = zeros (M, N);
VK = zeros (M, N);
P = zeros (M, N);
PSIK = zeros (M, N);
ZETAK = 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

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
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


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=UK;
end
vignesh is offline   Reply With Quote

Old   May 15, 2011, 10:56
Default lid driven Cavity using FVM
  #17
New Member
 
ainal
Join Date: May 2011
Posts: 2
Rep Power: 0
ainal is on a distinguished road
hello ,
I am working on a CFD problem and developing a Matlab code for "Lid driven cavity " problem using finite volume method and Symmetric Couple Gauss Seidel Scheme but i have some confusion in writing the code in Matlab so can any one help me in this ...

I shall be very grateful if anyone can share this
ainal is offline   Reply With Quote

Old   April 15, 2013, 00:15
Default
  #18
New Member
 
TAMILNADU
Join Date: Apr 2013
Posts: 2
Rep Power: 0
ananthulasriram is on a distinguished road
bro do u have any unsteady code for the same problem.. ??
ananthulasriram is offline   Reply With Quote

Old   April 14, 2014, 06:23
Default Lid Driven Cavity using Staggered grid
  #19
New Member
 
sankha
Join Date: Mar 2014
Posts: 1
Rep Power: 0
fireblade is on a distinguished road
Hi,
i have done LDC in collocated grid in stream function vorticity method...right now trying to do it in staggered grid in primitive variable method FDM....can somebody send me the code so that i can verify my results.. sm4fn@hotmail.com
i'm using SIMPLE
fireblade is offline   Reply With Quote

Old   July 21, 2014, 10:29
Default
  #20
New Member
 
Bob Man
Join Date: Jul 2014
Posts: 1
Rep Power: 0
Boobale is on a distinguished road
Quote:
Originally Posted by idris View Post
ok, here it is

clear all;clf;clc;
re=100; % for higher reynolds number(>500-1000), use under relaxation
rex=1; % underrelaxation, 0<rex<1
rexb=rex;

nx=51; dx=1/(nx-1); %100x100 grid (101*101 nodes)
ny=51; dy=1/(ny-1);

fac=2*(1/dx^2+1/dy^2);
Y=1;
uo=1;
nu=Y/(uo*re);


% initial condition

p(1:nx,1:ny)=0;
o(1:nx,1:ny)=0;
psi(1:nx,1:ny)=0;
omega(1:nx,1:ny)=0;
uact(1:nx,1:ny)=0;
vact(1:nx,1:ny)=0;
u(1:nx,1:ny)=0;
v(1:nx,1:ny)=0;
u(1:nx,ny)=1;

for iteration =1:2000%iteration

disp(['iteration= ',int2str(iteration)])
% stream function(internal nodes)

for i=2nx-1)
for j=2ny-1)
pold(i,j)=p(i,j);
f1=1/fac*(o(i,j)+(p(i+1,j)+p(i-1,j))/dx^2+(p(i,j+1) ...
+p(i,j-1))/dy^2);
p(i,j)=p(i,j)+rex*(f1-p(i,j));
end
end

% Vorticity (Boundary nodes)
for j=1:ny
fab=-2*p(2,j)/dx^2; %left side(AB)
o(1,j)=o(1,j)+rexb*(fab-o(1,j));
fcd=-2*p(nx-1,j)/dx^2; %right side(CD)
o(nx,j)=o(nx,j)+rexb*(fcd-o(nx,j));
end

for i=1:nx
fad=-2*p(i,2)/dy^2; %bottom side(AD)
o(i,1)=o(i,1)+rexb*(fad-o(i,1));
fbc=-(2*p(i,ny-1)+2*uo*dy)/dy^2; %top side (BC)
o(i,ny)=o(i,ny)+rexb*(fbc-o(i,ny));
end

%vorticity internal nodes

for i=2nx-1)
for j=2ny-1)
oold(i,j)=o(i,j);
f2=1/fac*((o(i+1,j)+o(i-1,j))/dx^2+(o(i,j+1)+o(i,j-1))/dy^2 ...
-re*(p(i,j+1)-p(i,j-1))*(o(i+1,j)-o(i-1,j))/(4*dx*dy)...
+re*(p(i+1,j)-p(i-1,j))*(o(i,j+1)-o(i,j-1))/(4*dx*dy));

o(i,j)=o(i,j)+rex*(f2-o(i,j));


end
end

% Calculate the velocity (u and v)

for i=2nx-1)
for j=2ny-1)
u(i,j)=(p(i,j+1)-p(i,j-1))/(2*dy);
v(i,j)=-(p(i+1,j)-p(i-1,j))/(2*dx);
end
end

% rearrange in cartesian coordinate
%(which originally in matrix coordinate)
for i=1:nx
for j=1:ny
k=i;
l=j;
psi(l,k)=p(i,j);
omega(l,k)=o(i,j);
uact(l,k)=u(i,j);
vact(l,k)=v(i,j);
end
end
figure(1)% Streamline plot with number
Z = psi(1:ny,1:nx);
X = linspace(0,1,size(Z,2));
Y = linspace(0,1,size(Z,1));
[c,h] = contour(X,Y,Z);

axis equal
axis([0 1 0 1])
drawnow
end

hope you can further developed these code by yourself. It's up to you what you are going to do about it :-)

Can anyone explain why the boundary conditions for the vorticity are what they are? I mean what kind of approximation is that for the derivative of psi?

Also what scheme is used for calculating the internal nodes of the vorticity?

Thanks alot for the help
Boobale 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
finite volume method - lid driven cavity ainal Main CFD Forum 0 May 15, 2011 11:06
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
Lid driven cavity flow KK FLUENT 1 December 16, 2009 19:10


All times are GMT -4. The time now is 11:27.