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

lid driven cavity matlab code

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

Reply
 
LinkBack Thread Tools Display Modes
Old   August 3, 2012, 15:43
Unhappy lid driven cavity matlab code
  #1
New Member
 
anna
Join Date: Aug 2012
Posts: 1
Rep Power: 0
anna_simons is on a distinguished road
hi, my name is anna, i am a CFD student
I developed a matlab code to solve 2D lid-driven cavity problem with finite volume scheme on staggered grid.
here is my code, but it don't work correctly, u loop don't converge for small residuals, and the answers aren't so good,
help me, tnx

clear all
%Problem Parameters
Re=100; %Reynolds No
Uw=1;
%Grid Properties
np=80; %grid No
dx=1/((np-2)/2);
dy=1/((np-2)/2);
%Initialization
Pstar=zeros(np,np);
Ustar=zeros(np,np);
Vstar=zeros(np,np);
convergence=zeros(1,2000);
k=0;
%under relaxation factor
relax=.7;
rel=.7;
Ustar(:,np-1)=Uw;
qsum=1;
%------------------------------------------------
%SIMPLE LOOP
while sqrt(qsum)>10^-2
k=k+1;
%display(k)
U=zeros(np,np);
V=zeros(np,np);
%Calculation of U field from Pstr and Ustr and Vsatr
U(:,np-1)=Uw;
residual1=1;
residual2=1;
while residual1>10^(-2)
if(residual1>residual2&&residual1>100)
display('U cal failed')
display(residual1)
break;
end
residual2=residual1;
for j=2.5: (np-3)/2)
for i=1.5: (np-3)/2)
b=dy*max(Ustar(2*i-1,2*j),0)+1/Re;
c=-((Vstar(2*i+2,2*j+1)+Vstar(2*i,2*j+1))/4)*dx+1/Re;
d=((Vstar(2*i+2,2*j-1)+Vstar(2*i,2*j-1))/4)*dx+1/Re;
e=(dy*max(-Ustar(2*i+3,2*j),0))/2+1/Re;
a=dy*abs(Ustar(2*i+1,2*j))-c-d+6/Re;
U(2*i+1,2*j)=(1-relax)*U(2*i+1,2*j)+relax*((b*U(2*i-1,2*j)+c*U(2*i+1,2*j+2)+d*U(2*i+1,2*j-2)+e*U(2*i+3,2*j)-dy*(Pstar(2*i+2,2*j)-Pstar(2*i,2*j)))/a);
end
end
delF=ones(1,((np)/2-3)*((np)/2-2));
s=1;
for j=2.5: (np-3)/2)
for i=1.5: (np-3)/2)
b=dy*max(Ustar(2*i-1,2*j),0)+1/Re;
c=-(Vstar(2*i+2,2*j+1)+Vstar(2*i,2*j+1))/4*dx+1/Re;
d=(Vstar(2*i+2,2*j-1)+Vstar(2*i,2*j-1))/4*dx+1/Re;
e=dy/2*max(-Ustar(2*i+3,2*j),0)+1/Re;
a=dy*abs(Ustar(2*i+1,2*j))+6/Re-c-d;
delF(s)=U(2*i+1,2*j)-((b*Ustar(2*i-1,2*j)+c*Ustar(2*i+1,2*j+2)+d*Ustar(2*i+1,2*j-2)+e*Ustar(2*i+3,2*j)-dy*(Pstar(2*i+2,2*j)-Pstar(2*i,2*j)))/a);
s=s+1;
end
end
sum=0;
for s=1:length(delF)
sum=sum+delF(s)^2;
end
residual1=sqrt(sum/length(delF));
display(residual1)
end
%Calculation of V field from Pstr and Ustr and Vsatr
while residual1>10^(-2)
if(residual1>residual2&&residual1>100)
display('V cal failed')
display(residual1)
break;
end
residual2=residual1;
for j=1.5: ((np-1)-2)/2)
for i=2.5: ((np-1)-2)/2)
B=dx*max(Vstar(2*i,2*j-1),0)+1/Re;
C=-((Ustar(2*i+1,2*j+2)+Ustar(2*i+1,2*j))/4)*dy+1/Re;
D=((Ustar(2*i-1,2*j+2)+Ustar(2*i-1,2*j))/4)*dy+1/Re;
E=(dx*max(-Vstar(2*i,2*j+3),0))/2+1/Re;
A=dx*abs(Vstar(2*i,2*j+1))+6/Re-C-D;
V(2*i,2*j+1)=(1-relax)* V(2*i,2*j+1)+relax*((B*V(2*i,2*j-1)+C*V(2*i+2,2*j+1)+D*V(2*i-2,2*j+1)+E*V(2*i,2*j+3)-dx*(Pstar(2*i,2*j+2)-Pstar(2*i,2*j)))/A);
end
end
delF=ones(1,((np)/2-3)*((np)/2-2));
s=1;
for j=1.5: ((np-1)-2)/2)
for i=2.5: ((np-1)-2)/2)
B=dx*max(Vstar(2*i,2*j-1),0)+1/Re;
C=-(Ustar(2*i+1,2*j+2)+Ustar(2*i+1,2*j))/4*dy+1/Re;
D=(Ustar(2*i-1,2*j+2)+Ustar(2*i-1,2*j))/4*dy+1/Re;
E=dx/2*max(-Vstar(2*i,2*j+3),0)+1/Re;
A=dx*abs(Vstar(2*i,2*j+1))+6/Re-C-D;
delF(s)= V(2*i,2*j+1)-((B*V(2*i,2*j-1)+C*V(2*i+2,2*j+1)+D*V(2*i-2,2*j+1)+E*V(2*i,2*j+3)-dx*(Pstar(2*i,2*j+2)-Pstar(2*i,2*j)))/A);
s=s+1;
end
end
sum=0;
for s=1:length(delF)
sum=sum+delF(s)^2;
end
residual1=sqrt(sum/length(delF));
end
Ustar=U;
Vstar=V;
Pcor=zeros(np,np);
%Calculation of P correction from Pstr and Ustr and Vstr
residual1=1;
residual2=1;
while residual1>10^(-2)
if(residual1>residual2&&residual1>100)
display('Pcor cal failed')
break;
end
residual2=residual1;
for j=2.5: ((np-1)-2)/2)
for i=2.5: ((np-1)-2)/2)
c2=dy/(dy*abs(Ustar(2*i+1,2*j))+4/Re+(Vstar(2*i+2,2*j+1)+Vstar(2*i,2*j+1)-Vstar(2*i+2,2*j-1)-Vstar(2*i,2*j-1))*dx/4);
c3=dy/(dy*abs(Ustar(2*i-1,2*j))+4/Re+(Vstar(2*i-2,2*j+1)+Vstar(2*i,2*j+1)-Vstar(2*i-2,2*j-1)-Vstar(2*i,2*j-1))*dx/4);
c4=dx/(dx*abs(Vstar(2*i,2*j+1))+4/Re+(Ustar(2*i+1,2*j+2)+Ustar(2*i+1,2*j)-Ustar(2*i-1,2*j+2)-Ustar(2*i-1,2*j))*dy/4);
c5=dx/(dx*abs(Vstar(2*i,2*j-1))+4/Re+(Ustar(2*i+1,2*j-2)+Ustar(2*i+1,2*j)-Ustar(2*i-1,2*j-2)-Ustar(2*i-1,2*j))*dy/4);
c6=-(U(2*i+1,2*j)-U(2*i-1,2*j)+V(2*i,2*j+1)-V(2*i,2*j-1));
c1=c2+c3+c4+c5;
if(i==2.5&&j==2.5)
Pcor(2*i,2*j)=(1-rel)* Pcor(2*i,2*j)+rel*((c2*Pcor(2*i+2,2*j)+c4*Pcor(2*i ,2*j+2)+c6)/c1);
elseif(i==2.5&&j==(((np-1)-2)/2))
Pcor(2*i,2*j)=(1-rel)* Pcor(2*i,2*j)+rel*((c2*Pcor(2*i+2,2*j)+c5*Pcor(2*i ,2*j-2)+c6)/c1);
elseif(i==(((np-1)-2)/2)&&j==2.5)
Pcor(2*i,2*j)=(1-rel)* Pcor(2*i,2*j)+rel*((c3*Pcor(2*i-2,2*j)+c4*Pcor(2*i,2*j+2)+c6)/c1);
elseif(i==(((np-1)-2)/2)&&j==(((np-1)-2)/2))
Pcor(2*i,2*j)=(1-rel)* Pcor(2*i,2*j)+rel*((c3*Pcor(2*i-2,2*j)+c4*Pcor(2*i,2*j+2)+c6)/c1);
elseif(i==2.5)
Pcor(2*i,2*j)=(1-rel)* Pcor(2*i,2*j)+rel*((c2*Pcor(2*i+2,2*j)+c4*Pcor(2*i ,2*j+2)+c5*Pcor(2*i,2*j-2)+c6)/c1);
elseif(i==(((np-1)-2)/2))
Pcor(2*i,2*j)=(1-rel)* Pcor(2*i,2*j)+rel*((c3*Pcor(2*i-2,2*j)+c4*Pcor(2*i,2*j+2)+c5*Pcor(2*i,2*j-2)+c6)/c1);
elseif(j==2.5)
Pcor(2*i,2*j)=(1-rel)* Pcor(2*i,2*j)+rel*((c2*Pcor(2*i+2,2*j)+c3*Pcor(2*i-2,2*j)+c4*Pcor(2*i,2*j+2)+c6)/c1);
elseif(j==(((np-1)-2)/2))
Pcor(2*i,2*j)=(1-rel)* Pcor(2*i,2*j)+rel*((c2*Pcor(2*i+2,2*j)+c3*Pcor(2*i-2,2*j)+c5*Pcor(2*i,2*j-2)+c6)/c1);
else
Pcor(2*i,2*j)=(1-rel)* Pcor(2*i,2*j)+rel*((c2*Pcor(2*i+2,2*j)+c3*Pcor(2*i-2,2*j)+c4*Pcor(2*i,2*j+2)+c5*Pcor(2*i,2*j-2)+c6)/c1);
end
end
end
delF=ones(1,((np)/2-5)*((np)/2-5));
s=1;
for j=2.5: ((np-1)-2)/2)
for i=2.5: ((np-1)-2)/2)
c2=dy/(dy*abs(Ustar(2*i+1,2*j))+4/Re+(Vstar(2*i+2,2*j+1)+Vstar(2*i,2*j+1)-Vstar(2*i+2,2*j-1)-Vstar(2*i,2*j-1))*dx/4);
c3=dy/(dy*abs(Ustar(2*i-1,2*j))+4/Re+(Vstar(2*i-2,2*j+1)+Vstar(2*i,2*j+1)-Vstar(2*i-2,2*j-1)-Vstar(2*i,2*j-1))*dx/4);
c4=dx/(dx*abs(Vstar(2*i,2*j+1))+4/Re+(Ustar(2*i+1,2*j+2)+Ustar(2*i+1,2*j)-Ustar(2*i-1,2*j+2)-Ustar(2*i-1,2*j))*dy/4);
c5=dx/(dx*abs(Vstar(2*i,2*j-1))+4/Re+(Ustar(2*i+1,2*j-2)+Ustar(2*i+1,2*j)-Ustar(2*i-1,2*j-2)-Ustar(2*i-1,2*j))*dy/4);
c6=-(U(2*i+1,2*j)-U(2*i-1,2*j)+V(2*i,2*j+1)-V(2*i,2*j-1));
c1=c2+c3+c4+c5;
if(i==2.5&&j==2.5)
delF(s)=Pcor(2*i,2*j)-((c2*Pcor(2*i+2,2*j)+c4*Pcor(2*i,2*j+2)+c6)/c1);
elseif(i==2.5&&j==(((np-1)-2)/2))
delF(s)=Pcor(2*i,2*j)-((c2*Pcor(2*i+2,2*j)+c5*Pcor(2*i,2*j-2)+c6)/c1);
elseif(i==(((np-1)-2)/2)&&j==2.5)
delF(s)=Pcor(2*i,2*j)-((c3*Pcor(2*i-2,2*j)+c4*Pcor(2*i,2*j+2)+c6)/c1);
elseif(i==(((np-1)-2)/2)&&j==(((np-1)-2)/2))
delF(s)=Pcor(2*i,2*j)-((c3*Pcor(2*i-2,2*j)+c4*Pcor(2*i,2*j+2)+c6)/c1);
elseif(i==2.5)
delF(s)=Pcor(2*i,2*j)-((c2*Pcor(2*i+2,2*j)+c4*Pcor(2*i,2*j+2)+c5*Pcor(2* i,2*j-2)+c6)/c1);
elseif(i==(((np-1)-2)/2))
delF(s)=Pcor(2*i,2*j)-((c3*Pcor(2*i-2,2*j)+c4*Pcor(2*i,2*j+2)+c5*Pcor(2*i,2*j-2)+c6)/c1);
elseif(j==2.5)
delF(s)=Pcor(2*i,2*j)-((c2*Pcor(2*i+2,2*j)+c3*Pcor(2*i-2,2*j)+c4*Pcor(2*i,2*j+2)+c6)/c1);
elseif(j==(((np-1)-2)/2))
delF(s)=Pcor(2*i,2*j)-((c2*Pcor(2*i+2,2*j)+c3*Pcor(2*i-2,2*j)+c5*Pcor(2*i,2*j-2)+c6)/c1);
else
delF(s)=Pcor(2*i,2*j)-((c2*Pcor(2*i+2,2*j)+c3*Pcor(2*i-2,2*j)+c4*Pcor(2*i,2*j+2)+c5*Pcor(2*i,2*j-2)+c6)/c1);
end
s=s+1;
end
end
sum=0;
for s=1:length(delF)
sum=sum+delF(s)^2;
end
residual1=sqrt(sum/length(delF));
end
Pcor(3,: ) = Pcor(5,;
Pcor(np-1, = Pcor(np-3, : );
Pcor(:,np-1)=Pcor(:,np-3);
Pcor(:,3)=Pcor(:,5);
Ucor=zeros(np,np);
% calculate Ucor from Pcor
for j=2.5: ((np-1)-2)/2)
for i=1.5: ((np-1)-2)/2)
c=-(V(2*i+2,2*j+1)+V(2*i,2*j+1))/4*dx+1/Re;
d=(V(2*i+2,2*j-1)+V(2*i,2*j-1))/4*dx+1/Re;
a=dy*abs(U(2*i+1,2*j))+6/Re-c-d;
Ucor(2*i+1,2*j)=-(dy/a)*(Pcor(2*i+2,2*j)-Pcor(2*i,2*j));
end
end
% calculate Vcor from Pcor
Vcor=zeros(np,np);
for j=1.5: ((np-3)/2)
for i=2.5: ((np-3)/2)
C=-(U(2*i+1,2*j+2)+U(2*i+1,2*j))/4*dy+1/Re;
D=(Ustar(2*i-1,2*j+2)+U(2*i-1,2*j))/4*dy+1/Re;
A=dx*abs(V(2*i,2*j+1))+6/Re-C-D;
Vcor(2*i,2*j+1)=-(dx/A)*(Pcor(2*i,2*j+2)-Pcor(2*i,2*j));
end
end
U=Ustar+Ucor;
V=Vstar+Vcor;
Pstar=Pstar+.3*Pcor;
Q=zeros(np,np);
qsum=0;
for j=2.5: ((np-3)/2)
for i=2.5: ((np-3)/2)
Q(2*i,2*j)=Ustar(2*i+1,2*j)-Ustar(2*i-1,2*j)+Vstar(2*i,2*j+1)-Vstar(2*i,2*j-1);
qsum=qsum+Q(2*i,2*j)^2;
end
end
convergence(k)=sqrt(qsum);
display(qsum)
end

Last edited by anna_simons; August 3, 2012 at 16:00.
anna_simons is offline   Reply With Quote

Old   August 3, 2012, 16:21
Default
  #2
Senior Member
 
Join Date: Aug 2011
Posts: 251
Rep Power: 7
leflix is on a distinguished road
Hi Anna,

waouh!! I'm not sure you will find someone here to debugg your code but I may be wrong.... ;-)
What you could do , should be to find here one code written in matlab for the same problem (I have seen some over there or you can find some on internet).

here are few
Matlab code for pipe flow
http://math.mit.edu/cse/codes/mit18086_navierstokes.pdf
http://www-math.mit.edu/cse/codes/mi...navierstokes.m

So get one and try to understand where you have made a mistake.

good luck
leflix is offline   Reply With Quote

Old   August 3, 2012, 17:18
Default
  #3
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,660
Rep Power: 23
FMDenaro will become famous soon enough
I can not debug all lines of your code but I suggest you to run for 1 time-step and check after if the divergence-free constraint in ensured in all the cells.
Good luck
FMDenaro is online now   Reply With Quote

Old   August 8, 2012, 11:00
Default
  #4
New Member
 
Join Date: Jul 2012
Posts: 5
Rep Power: 5
Eifoehn4 is on a distinguished road
In your SOR-Iteration of u,v and p:

Do you take the "old" and the "new" value into account?

Regards
Eifoehn4 is offline   Reply With Quote

Reply

Tags
cavity problem

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
Boundary condition for 2d lid driven cavity using ghost cells quarkz Main CFD Forum 9 January 20, 2013 06:54
lid-driven cavity in matlab using BiCGStab Don456 Main CFD Forum 1 January 19, 2012 16:00
Lid Driven Cavity using Ghost Cell Method and in C++ illuminati5288 Main CFD Forum 0 August 12, 2011 22:05
Lack of Recirculation for Lid Driven Cavity Flow klw Main CFD Forum 8 May 21, 2011 04:57
help wanted abt lid driven cavity fazle rabbi ahad Main CFD Forum 2 February 20, 2007 09:23


All times are GMT -4. The time now is 09:17.