# Need help debugging lid driven cavity flow that I tried to develop

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

 October 17, 2017, 13:00 Need help debugging lid driven cavity flow that I tried to develop #1 New Member   Antony Join Date: Oct 2017 Posts: 2 Rep Power: 0 Hi all, So here is the story, I had a 2D SIMPLE ( semi-implicit method for pressure-linked equations) code for solving incompressible navier stokes equation in Matlab. My job is to develop a two phase flow model, so I started by adding density into my equations. but the problem is, now that I added the density into my code, it is not giving me the right answer for lid driven cavity problem with constant density and I have no idea why. I tried several boundary conditions, I tried changing lots of thing (which might have added more bugs in my code). Could you take a quick look at it and let me know if you can find any bug? Thank you A. here is the code: working_cavity.zip P.S. If the bug is found, this is a good code for others to use Also, sorry for bad coding

 October 17, 2017, 13:10 #2 New Member   Antony Join Date: Oct 2017 Posts: 2 Rep Power: 0 Here is the code: Again sorry that it has few comments: function working_cavity nx=20; ny=20; x_max=1; x_min=0.0; y_max=1.0; y_min=0.0; Urfu=.7; %under relaxation u Urfv=.7; Urfp=.5; deltax = (x_max-x_min)/(nx-1); deltay = (y_max-y_min)/(ny-1); Viscos=1e-2; ro_fluid=1.0; Uin=1; Re=ro_fluid*Uin*(x_max-x_min)/Viscos Asu=deltax; Awe=deltay; Aea=Awe; Ano=Asu; Vol=Asu*Awe; Sormax =0.00001; NswpU=3; NswpV=2; NswpP=10; t_final=0.2; deltat=0.1; Maxts=round(t_final/deltat+2); Great=10e30; fileID = fopen('results.dat','w'); fmt = '%f %f %f %f %f\n'; fprintf(fileID, 'VARIABLES = "x","y","P","U","V" \n'); for i = 1:nx-1 %changed to -1 for j = 1:ny-1 if (i~=nx-1 && j~=ny-1) XC(i,j) = x_min + (i-0.5)*deltax; YC(i,j) = y_min + (j-0.5)*deltay; end X(i,j)=x_min + (i-1)*deltax; Y(i,j)=y_min + (j-1)*deltay; end end U_fluid=zeros(nx,ny); V_fluid=zeros(nx,ny); U_fluid_old=zeros(nx,ny); V_fluid_old=zeros(nx,ny); Pp=zeros(nx,ny); P=zeros(nx,ny); Du=zeros(nx,ny); Dv=zeros(nx,ny); Su=zeros(nx,ny); Sp=zeros(nx,ny); for I=2:nx-1 U_fluid(I,ny-1)=Uin; U_fluid(I,ny)=Uin; end max_iter=100; for N=2:Maxts-1 t=(N-1)*deltat; Kr=0; Source1=10.0; for I=1:nx for J=1:ny index1 = I + nx*(J-1); Densit(I,J)=ro_fluid; Densit_old(I,J)=ro_fluid; end end while (Source1 > Sormax) Kr=Kr+1; for I=1:nx for J=1:ny An(I,J)=0.0; Ap(I,J)=0.0; As(I,J)=0.0; Ae(I,J)=0.0; Aw(I,J)=0.0; Su(I,J)=0.0; Sp(I,J)=0.0; Du(I,J)=Great; Dv(I,J)=Great; end end % C ------------------Calc U --------------------- for I=3:nx-1 for J=2:ny-1 % C -----Calculate Convection Coefficients-------------- Ce=Densit(I,J)*.5*(U_fluid(I,J)+U_fluid(I+1,J))*Ae a; % Check density Cw=Densit(I-1,J)*.5*(U_fluid(I,J)+U_fluid(I-1,J))*Awe; % Check density Cs=0.25*(Densit(I,J)+Densit(I,J-1)+Densit(I-1,J)+Densit(I-1,J-1))*.5*(V_fluid(I,J)+V_fluid(I-1,J))*Asu; Cn=0.25*(Densit(I,J)+Densit(I,J+1)+Densit(I-1,J)+Densit(I-1,J+1))*.5*(V_fluid(I,J+1)+V_fluid(I-1,J+1))*Ano; % C -----Calculate Diffusion Coefficients Ds=Viscos*Asu/(deltay); De=Viscos*(Aea/(deltax)); Dw=Viscos*(Awe/(deltax)); Dn=Viscos*Ano/(deltay); % C -----Assemble Main Coefficients An(I,J)= Dn+max([-Cn/2.0, 0.0]); As(I,J)= Ds+max([ Cs/2.0, 0.0]); Ae(I,J)= De+max([-Ce/2.0, 0.0]); Aw(I,J)= Dw+max([ Cw/2.0, 0.0]); % C -----Calculate Coefficients Of Source Terms Gue=(U_fluid(I,J)-U_fluid(I+1,J))/(-1*deltax); Guw=(U_fluid(I,J)-U_fluid(I-1,J))/(deltax); Gvn=(V_fluid(I,J+1)-V_fluid(I-1,J+1))/(deltax); Gvs=(V_fluid(I,J)-V_fluid(I-1,J))/(deltax); Spp=-Awe*(P(I,J)-P(I-1,J)); Gen=0.0; Gen=(Viscos*Gue-Viscos*Guw)*Awe+(Viscos*Gvn-Viscos*Gvs)*Ano; Smp=Cn-Cs+Ce-Cw; Cp=max(0.0,Smp); Sp(I,J)=-Cp; Su(I,J)=Cp*U_fluid(I,J)+Spp; Su(I,J)=Su(I,J)+Gen; end end % C %%%%%%%%%% BOUNDARY U %%%%%%%%% (learn more) % C TOP WALL J=ny-1; for I=3:nx-1 U_fluid(I,J)=Uin; U_fluid(I,J+1)=U_fluid(I,J); Tmult=Viscos*deltax/deltay; %% find if it is right try Viscos*deltax/deltay Sp(I,J)=Sp(I,J)-Tmult; Su(I,J)=Su(I,J)+Tmult*U_fluid(I,J+1); An(I,J)=0.0; end % C SOUTH WALL J=2; for I=3:nx-1 U_fluid(I,J)=0.0; U_fluid(I,J-1)=U_fluid(I,J); Tmult=Viscos*deltax/deltay; Sp(I,J)=Sp(I,J)-Tmult; Su(I,J)=Su(I,J)+Tmult*U_fluid(I,J-1); As(I,J)=0.0; end % C East WALL I=nx-1;%play for J=2:ny-1 U_fluid(nx,J)=0.0; Aw(I,J)=0.0; As(I,J)= 0.0; Ae(I,J)= 0.0; An(I,J)=0.0; % Su(I,J)=0.0; % Sp(I,J)=-Great; end % C West WALL I=3; for J=2:ny-1 U_fluid(2,J)=0.0; Aw(I,J)=0.0; As(I,J)= 0.0; Ae(I,J)= 0.0; An(I,J)=0.0; % Su(I,J)=0.0; % Sp(I,J)=-Great; end Resoru=0.0; for I=3:nx-1 %double checked for J=2:ny-1 Ap(I,J)=An(I,J)+As(I,J)+Ae(I,J)+Aw(I,J)-Sp(I,J); % C ----- transient Su(I,J)=Su(I,J)+0.5*(Densit_old(I-1,J)+Densit_old(I,J))*Vol/deltat*U_fluid_old(I,J); Ap(I,J)=Ap(I,J)+0.5*(Densit(I-1,J)+Densit(I,J))*Vol/deltat; Du(I,J)=Ap(I,J); Resor=An(I,J)*U_fluid(I,J+1)+As(I,J)*U_fluid(I,J-1)+Ae(I,J)*U_fluid(I+1,J)+Aw(I,J)*U_fluid(I-1,J)-Ap(I,J)*U_fluid(I,J)+Su(I,J); Resoru=Resoru+abs(Resor); % C --- Under-Relaxation Ap(I,J)=Ap(I,J)/Urfu; Su(I,J)=Su(I,J)+(1.0-Urfu)*Ap(I,J)*U_fluid(I,J); Du(I,J)=Du(I,J)*Urfu; %Simple % Du(I,J)=Du(I,J)*(1.0-Urfu)/Urfu-Sp(I,J); %!Simplec end end for Nu=1:NswpU U_fluid=Liner_Solver(3,2,U_fluid,nx,ny,An,As,Ae,Aw ,Su,Ap); end % C ------------------ End Calc U --------------------- for I=1:nx for J=1:ny An(I,J)=0.0; Ap(I,J)=0.0; As(I,J)=0.0; Ae(I,J)=0.0; Aw(I,J)=0.0; Su(I,J)=0.0; Sp(I,J)=0.0; end end % C ------------------ Calc V --------------------- for I=2:nx-1 for J=3:ny-1 Ce=0.25*(Densit(I,J)+Densit(I+1,J)+Densit(I,J-1)+Densit(I+1,J-1))*0.5*(U_fluid(I+1,J)+U_fluid(I+1,J-1))*Aea; Cw=0.25*(Densit(I,J)+Densit(I-1,J)+Densit(I,J-1)+Densit(I-1,J-1))*0.5*(U_fluid(I,J)+U_fluid(I,J-1))*Awe; Cn=Densit(I,J)*0.5*(V_fluid(I,J)+V_fluid(I,J+1))*A no; Cs=Densit(I,J-1)*0.5*(V_fluid(I,J)+V_fluid(I,J-1))*Asu; Dn=Viscos*Ano/deltay; Ds=Viscos*Asu/deltay; De=Viscos*Aea/deltax; Dw=Viscos*Awe/deltax; Spp=-Ano*(P(I,J)-P(I,J-1));%!+Gravity*.5*(Densit(I,J)+Densit(I,J-1))*Vol An(I,J)= Dn+max([-Cn/2.0, 0.0]); As(I,J)= Ds+max( [Cs/2.0, 0.0]); Ae(I,J)= De+max([-Ce/2.0, 0.0]); Aw(I,J)= Dw+max( [Cw/2.0, 0.0]); Smp=Cn-Cs+Ce-Cw; Gue=(U_fluid(I+1,J)-U_fluid(I+1,J-1))/(deltay); Guw=(U_fluid(I,J)-U_fluid(I,J-1))/(deltay); Gvn=(V_fluid(I,J)-V_fluid(I,J+1))/(-1*deltay); Gvs=(V_fluid(I,J)-V_fluid(I,J-1))/(deltay); Gen=0.0;%-Gravity*.5*(Densit(I,J)+Densit(I,J-1))*Vol; % Gen=Gen+(Viscos*Gue-Viscos*Guw)*Awe+(Viscos*Gvn-Viscos*Gvs)*Ano; Cp=max(0.0,Smp); Sp(I,J)=-Cp; Su(I,J)=Cp*V_fluid(I,J)+Spp; Su(I,J)=Su(I,J)+Gen; end end % C %%%%%%%%%% BOUNDARY V %%%%%%%%% % !-----North J=ny-1; for I=2:nx-1 V_fluid(I,ny)=0.0; Aw(I,J)=0.0; As(I,J)= 0.0; Ae(I,J)= 0.0; An(I,J)=0.0; % Su(I,J)=0.0; % Sp(I,J)=-Great; end % !-----South J=3; for I=2:nx-1 V_fluid(I,2)=0.0; Aw(I,J)=0.0; As(I,J)= 0.0; Ae(I,J)= 0.0; An(I,J)=0.0; % Su(I,J)=0.0; % Sp(I,J)=-Great; end % !-----West I=2; for J=3:ny-1 V_fluid(I,J)=0.0; V_fluid(I-1,J)=V_fluid(I,J); Tmult=Viscos*deltax/deltay; Sp(I,J)=Sp(I,J)-Tmult; Su(I,J)=Su(I,J)+Tmult*V_fluid(I-1,J); Aw(I,J)=0.0; end % !-----Eest I=nx-1; for J=3:ny-1 V_fluid(I,J)=0.0; V_fluid(I+1,J)=V_fluid(I,J); Tmult=Viscos*deltax/deltay; Sp(I,J)=Sp(I,J)-Tmult; Su(I,J)=Su(I,J)+Tmult*V_fluid(I+1,J); Ae(I,J)=0.0; end Resorv=0.0; for I=2:nx-1 for J=3:ny-1 Ap(I,J)=An(I,J)+As(I,J)+Ae(I,J)+Aw(I,J)-Sp(I,J); % C ! Transient Su(I,J)=Su(I,J)+0.5*(Densit_old(I,J-1)+Densit_old(I,J))*Vol/deltat*V_fluid_old(I,J); Ap(I,J)=Ap(I,J)+0.5*(Densit(I,J-1)+Densit(I,J))*Vol/deltat; Dv(I,J)=Ap(I,J); Resor=An(I,J)*V_fluid(I,J+1)+As(I,J)*V_fluid(I,J-1)+Ae(I,J)*V_fluid(I+1,J)+Aw(I,J)*V_fluid(I-1,J)-Ap(I,J)*V_fluid(I,J)+Su(I,J); Resorv=Resorv+abs(Resor); Ap(I,J)=Ap(I,J)/Urfv; Su(I,J)=Su(I,J)+(1.0-Urfv)*Ap(I,J)*V_fluid(I,J); Dv(I,J)=Dv(I,J)*Urfv; %!Simple % Dv(I,J)=Dv(I,J)*(1.0-Urfv)/Urfv-Sp(I,J);% !Simplec end end for Nv=1:NswpV V_fluid=Liner_Solver(2,3,V_fluid,nx,ny,An,As,Ae,Aw ,Su,Ap); end % C ------------------ Endof Calc V --------------------- for I=1:nx for J=1:ny An(I,J)=0.0; Ap(I,J)=0.0; As(I,J)=0.0; Ae(I,J)=0.0; Aw(I,J)=0.0; Su(I,J)=0.0; Sp(I,J)=0.0; end end % C ------------------ Calc P --------------------- (learn more) Resorm=0.0; for I=2:nx-1 for J=2:ny-1 % !----- North An(I,J)=-0.5*(Densit(I,J)+Densit(I,J+1))*(Ano)^2.0/Dv(I,J+1); G2starn=0.5*(Densit(I,J)+Densit(I,J+1))*V_fluid(I, J+1)*Ano; As(I,J)=-0.5*(Densit(I,J)+Densit(I,J-1))*(Asu)^2.0/Dv(I,J); G2stars=0.5*(Densit(I,J)+Densit(I,J-1))*V_fluid(I,J)*Asu; Ae(I,J)=-0.5*(Densit(I,J)+Densit(I+1,J))*(Awe)^2.0/Du(I+1,J); G1stare=0.5*(Densit(I,J)+Densit(I+1,J))*U_fluid(I+ 1,J)*Aea; Aw(I,J)=-0.5*(Densit(I,J)+Densit(I-1,J))*(Aea)^2.0/Du(I,J); G1starw=0.5*(Densit(I,J)+Densit(I-1,J))*U_fluid(I,J)*Aea; if J==ny-1 An(I,J)=0.0; G2starn=0.0; end % % !----- South if J==2 As(I,J)=0.0; G2stars=0.0; end % % !----- West if I==2 Aw(I,J)=0.0; G1starw=0.0; end % % !----- East if I==nx-1 Ae(I,J)=0.0; G1stare=0.0; end Smp=G1stare-G1starw+G2starn-G2stars; Su(I,J)=Smp+(Densit(I,J)-Densit_old(I,J))/deltat*Vol; Sp(I,J)=0.0; Resorm=Resorm+abs(Smp); end end for I=2:nx-1 for J=2:ny-1 Ap(I,J)=Ae(I,J)+Aw(I,J)+An(I,J)+As(I,J);%-Sp(I,J); end end for Np=1:NswpP Pp=Liner_Solver(2,2,Pp,nx,ny,An,As,Ae,Aw,Su,Ap); end % C ------------------ Endof Calc P --------------------- % C ------------------ Correction --------------------- % C -------------------- U for I=3:nx-1 for J=2:ny-1 U_fluid(I,J)=U_fluid(I,J)-Awe/Du(I,J)*(Pp(I,J)-Pp(I-1,J)); end end % C -------------------- V for I=2:nx-1 for J=3:ny-1 V_fluid(I,J)=V_fluid(I,J)-Asu/Dv(I,J)*(Pp(I,J)-Pp(I,J-1)); end end % C --------------------- Pressures (what is pressure refrence?) Ppref=Pp(2,2); for I=2:nx-1 for J=2:ny-1 P(I,J)=P(I,J)+Urfp*(Pp(I,J)-Ppref); Pp(I,J)=0.0; end end % C ------------------ Endof Correction --------------------- plot_everything(X,Y,XC,YC,nx,ny,U_fluid,V_fluid,P) ; Source1=max([Resoru,Resorv,Resorm]) if ( Kr > max_iter) % ! 500 is max itteration per time Kr for i=1:nx-2 for j=1:ny-2 Ug(i,j)=(U_fluid(i+1,j+1)+U_fluid(i+2,j+1))/2.0; Vg(i,j)=(V_fluid(i+1,j+1)+V_fluid(i+1,j+2))/2.0; Pp(i,j)=P(i+1,j+1); end end if N == Maxts-1 title=['Zone I=',num2str(nx-2),' J=',num2str(ny-2),' t="',num2str(t) ,'"\n']; fprintf(fileID, title); for i=1:nx-2 for j=1:ny-2 fprintf(fileID,fmt, [XC(i,j) YC(i,j) Pp(i,j) Ug(i,j) Vg(i,j)]); end end figure(4); streamline(XC(:,1)',YC(1,,Ug,Vg,rand(50,1),rand( 50,1),[0.1,1000]); axis image; axis([0,1,0,1]); end break end end end fclose(fileID); end function Phi=Liner_Solver(Istart,Jstart,Phi,nx,ny,An,As,Ae, Aw,Su,Ap) Jstm1=Jstart-1; AA(Jstm1)=0; %%-----Commence W-E Sweep for I=Istart:nx-1 CC(Jstm1)=Phi(I,Jstm1); %%-----Commence S-N Traverse for J=Jstart:ny-1 %%-----Assemble Toma Coefficients AA(J)=An(I,J); BB(J)=As(I,J); CC(J)=Ae(I,J)*Phi(I+1,J)+Aw(I,J)*Phi(I-1,J)+Su(I,J); DD(J)=Ap(I,J); %%-----Calculate Coefficients Of Recorrence Formula Term=1./(DD(J)-BB(J)*AA(J-1)); AA(J)=AA(J)*Term; CC(J)=(BB(J)*CC(J-1)+CC(J))*Term; end %%-----Obtain New Phi's for Jj=Jstart:ny-1 J=ny+Jstm1-Jj; Phi(I,J)=AA(J)*Phi(I,J+1)+CC(J); end end end function plot_everything(X,Y,XC,YC,nx,ny,U_fluid,V_fluid,P) ro=zeros(nx,ny); for i=1:nx-2 for j=1:ny-2 Ug(i,j)=(U_fluid(i+1,j+1)+U_fluid(i+2,j+1))/2.0; Vg(i,j)=(V_fluid(i+1,j+1)+V_fluid(i+1,j+2))/2.0; Pp(i,j)=P(i+1,j+1); end end %%%%%%%%%%%%%%%%%%%%%%% CHANGED figure(2) % quiver(X,Y,Ug,Vg); quiver(XC,YC,Ug,Vg); xlabel('x'); ylabel('y'); title('velocity') figure (3) % [sx sy] = meshgrid(X(10), [Y(4),Y(10),Y(15)]); % % stream2(X,Y,Ug,Vg,sx,sy) contourf(XC,YC,Pp) colorbar pause (0.1) end

 Tags density, lid driven cavity, simple