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 Last edited by anna_simons; August 3, 2012 at 16:00.

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

 August 3, 2012, 17:18 #3 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 1,660 Rep Power: 23 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 New Member   Join Date: Jul 2012 Posts: 5 Rep Power: 5 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 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 On 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 06:54 Don456 Main CFD Forum 1 January 19, 2012 16: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 09:23

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