[Antony_R]

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]

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