[anna_simons]

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

[leflix]

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
http://www.cfd-online.com/Forums/mai...pipe-flow.html
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

[FMDenaro]

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

[Eifoehn4]

In your SOR-Iteration of u,v and p:

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

Regards