# lid driven cavity matlab code

 User Name Remember Me Password
 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: 272 Rep Power: 15 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,793 Rep Power: 71 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: 184 Rep Power: 13 In your SOR-Iteration of u,v and p: Do you take the "old" and the "new" value into account? Regards

 Tags cavity problem

 Thread Tools Search this Thread Search this Thread: Advanced Search Display Modes Linear Mode

 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 OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post quarkz Main CFD Forum 9 January 20, 2013 05:54 Don456 Main CFD Forum 1 January 19, 2012 15:00 illuminati5288 Main CFD Forum 0 August 12, 2011 22:05 klw Main CFD Forum 8 May 21, 2011 04:57 fazle rabbi ahad Main CFD Forum 2 February 20, 2007 08:23

All times are GMT -4. The time now is 10:57.

 Contact Us - CFD Online - Privacy Statement - Top