# lid driven cavity matlab code

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

 August 3, 2012, 15:43 lid driven cavity matlab code #1 New Member   anna Join Date: Aug 2012 Posts: 1 Rep Power: 0 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 Hesam_Ami likes this. Last edited by anna_simons; August 3, 2012 at 16:00.

 August 3, 2012, 16:21 #2 Senior Member   Join Date: Aug 2011 Posts: 270 Rep Power: 14 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

 August 3, 2012, 17:18 #3 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,395 Rep Power: 67 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

 August 8, 2012, 11:00 #4 Senior Member     - Join Date: Jul 2012 Location: Germany Posts: 175 Rep Power: 12 In your SOR-Iteration of u,v and p: Do you take the "old" and the "new" value into account? Regards

 Tags cavity problem