CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

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

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   October 17, 2017, 12:00
Default 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
Antony_R is on a distinguished road
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
Antony_R is offline   Reply With Quote

Old   October 17, 2017, 12:10
Default
  #2
New Member
 
Antony
Join Date: Oct 2017
Posts: 2
Rep Power: 0
Antony_R is on a distinguished road
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
Antony_R is offline   Reply With Quote

Reply

Tags
density, lid driven cavity, simple

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
lid driven cavity varying results yasmil OpenFOAM Running, Solving & CFD 2 October 6, 2016 21:42
Lid driven cavity problem thegodfather Main CFD Forum 1 December 2, 2012 04:37
Lid Driven Cavity using Ghost Cell Method and in C++ illuminati5288 Main CFD Forum 0 August 12, 2011 22:05
code for SIMPLE algorithm - 2D Lid driven cavity flow problem - Collocated grid h_amooie OpenFOAM Programming & Development 0 June 9, 2011 11:17
[OpenFOAM] Paraview - Lid Driven Cavity kieranhood ParaView 0 February 13, 2010 16:28


All times are GMT -4. The time now is 13:53.