[houkensjtu]

hey guys
i have made a new matlab code on laminar flow.It uses a staggered grid and SIP for iteration method.If u r also a matlab user working on cfd,u can download it and just run it on your pc.
And also i hope u can make some advice.
Thank u

%%%%%%%%%%%%%%%%%2D-cartesian coordinate SIMPLE Method%%%%%%%%%%%%%%%%%%%%%

clc
clear all;
format long;
%initialization
errdata=zeros(1,100000);
nx=120;
ny=30;
%velcity field
u=zeros(ny+2,nx+1,3);
v=zeros(ny+1,nx+2,3);
ua=zeros(ny+2,nx+1);
va=zeros(ny+1,nx+2);
ua(16:17,40:42)=10^25;
va(15:17,41:42)=10^25;
%pressure
p=zeros(ny,nx,3);
%co. for momentum equation
aeu=zeros(ny+2,nx+1);
awu=zeros(ny+2,nx+1);
asu=zeros(ny+2,nx+1);
anu=zeros(ny+2,nx+1);
apu=zeros(ny+2,nx+1);
aev=zeros(ny+1,nx+2);
awv=zeros(ny+1,nx+2);
asv=zeros(ny+1,nx+2);
anv=zeros(ny+1,nx+2);
apv=zeros(ny+1,nx+2);
%co. for SIP method
ls=zeros(ny+2,nx+1);
lw=zeros(ny+2,nx+1);
un=zeros(ny+2,nx+1);
ue=zeros(ny+2,nx+1);
lpr=zeros(ny+2,nx+1);
res=zeros(ny+2,nx+1);
resl=0;
%co. for pressure correction equation
aepc=zeros(ny,nx);
awpc=zeros(ny,nx);
anpc=zeros(ny,nx);
aspc=zeros(ny,nx);
appc=zeros(ny,nx);
pc=zeros(ny+2,nx+2,3);
%mass unbalance
mp=zeros(ny,nx);
%dimensionless density,grid size,viscosity
raw=1;
dx=1/ny;
dy=1/ny;
d=0.1;
alpha=0.0002;
%the inlet velocity=RE
for i=2:ny+1
u(i,:,=-6*((0.5/ny+(i-2)/ny)-0.5)^2+1.5;
end
u(16:17,40:42)=0;
v(15:17,41:42)=0;
for j=1:nx
p(:,j,=(1.2*nx/ny)/nx*(1-j);
end
%at this time mp is all zero
for i=1:ny
for j=1:nx
mp(i,j)=raw*dx*((u(i,j+1,1)-u(i,j,1))+(v(i+1,j,1)-v(i,j,1)));
end
end
%main loop
for q=1:100000
%%%%%%%%%%%%%%%%%%%%%%%%%momentum equation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ll=1:1
u(16:17,40:42)=0;
v(15:17,41:42)=0;
for i=2:ny+1
for j=2:nx
if i==2
aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
asu(i,j)=0;%d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
%awu(:,2)=0;
%aeu(:,nx)=0;
apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j)+0.6*d x/u(i,j,1);%(dy*(1.2*nx/ny)/nx+0.1*(u(i+1,j,1)-u(i,j,1)))/u(i,j,1);
%u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i-1,j-1,1)-p(i-1,j,1))*dy+ua(i,j)*u(i,j,1))/(apu(i,j)+ua(i,j));
else if i==ny+1
aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
anu(i,j)=0;%d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
%awu(:,2)=0;
%aeu(:,nx)=0;
apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j)+0.6*d x/u(i,j,1);%+(dy*(1.2*nx/ny)/nx+0.1*(u(i-1,j,1)-u(i,j,1)))/u(i,j,1);
%u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i-1,j-1,1)-p(i-1,j,1))*dy+ua(i,j)*u(i,j,1))/(apu(i,j)+ua(i,j));
else
aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
%awu(:,2)=0;
%aeu(:,nx)=0;
apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j);
%u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i-1,j-1,1)-p(i-1,j,1))*dy+ua(i,j)*u(i,j,1))/(apu(i,j)+ua(i,j));
end
end
end
end
%SIP METHOD for u velocity
su=zeros(ny+2,nx+1);
urf=0.02;
%co. for SIP method
ls=zeros(ny+2,nx+1);
lw=zeros(ny+2,nx+1);
un=zeros(ny+2,nx+1);
ue=zeros(ny+2,nx+1);
lpr=zeros(ny+2,nx+1);
res=zeros(ny+2,nx+1);
resl=0;
for i=2:ny+1
for j=2:nx
apu(i,j)=-(apu(i,j)+ua(i,j))/0.8;
su(i,j)=0.2*apu(i,j)*u(i,j,1)+dy*(p(i-1,j,1)-p(i-1,j-1,1))-ua(i,j)*u(i,j,1);
end
end
for j=2:nx
for i=2:ny+1
lw(i,j)=awu(i,j)/(1+0.9*un(i,j-1));
ls(i,j)=asu(i,j)/(1+0.9*ue(i-1,j));
p1=0.9*lw(i,j)*un(i,j-1);
p2=0.9*ls(i,j)*ue(i-1,j);
lpr(i,j)=1/(apu(i,j)+p1+p2-lw(i,j)*ue(i,j-1)-ls(i,j)*un(i-1,j)+10^-25);
un(i,j)=(anu(i,j)-p1)*lpr(i,j);
ue(i,j)=(aeu(i,j)-p2)*lpr(i,j);
end
end
for l=1:3
resl=0;
for i=2:ny+1
for j=2:nx
res(i,j)=su(i,j)-(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+apu(i,j)*u(i,j,1));
resl=resl+abs(res(i,j));
res(i,j)=(res(i,j)-ls(i,j)*res(i-1,j)-lw(i,j)*res(i,j-1))*lpr(i,j);
end
end
if l==1
resor=resl;
end
% rsm=resl/(resor+10^-25)
resl;
for i=(ny+1):-1:2
for j=nx:-1:2
res(i,j)=res(i,j)-un(i,j)*res(i+1,j)-ue(i,j)*res(i,j+1);
u(i,j,1)=u(i,j,1)+res(i,j);
end
end
end
%u(:,nx,1)=u(:,nx,1)*(sum(u(:,1,1))/sum(u(:,nx,1)));
u(:,nx+1,1)=u(:,nx,1);%for a outlet condition
for i=2:ny
for j=2:nx+1
if j==2
aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
awv(i,j)=0;%2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j-1,1)+0.25*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
else if j==nx+1
aev(i,j)=0;%2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j,1)+0.25*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
else
aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
end
end
v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+(p(i-1,j-1,1)-p(i,j-1,1))*dy+va(i,j)*v(i,j,1))/(apv(i,j)+va(i,j));
end
end
%SIP METHOD for v velocity
sv=zeros(ny+1,nx+2);
urf=0.02;
%co. for SIP method
ls=zeros(ny+1,nx+2);
lw=zeros(ny+1,nx+2);
un=zeros(ny+1,nx+2);
ue=zeros(ny+1,nx+2);
lpr=zeros(ny+1,nx+2);
res=zeros(ny+1,nx+2);
resl=0;
for i=2:ny
for j=2:nx+1
apv(i,j)=-(apv(i,j)+va(i,j))/0.8;
sv(i,j)=0.2*apv(i,j)*v(i,j,1)+dx*(p(i,j-1,1)-p(i-1,j-1,1))-va(i,j)*v(i,j,1);
end
end
for j=2:nx+1
for i=2:ny
lw(i,j)=awv(i,j)/(1+0.9*un(i,j-1));
ls(i,j)=asv(i,j)/(1+0.9*ue(i-1,j));
p1=0.9*lw(i,j)*un(i,j-1);
p2=0.9*ls(i,j)*ue(i-1,j);
lpr(i,j)=1/(apv(i,j)+p1+p2-lw(i,j)*ue(i,j-1)-ls(i,j)*un(i-1,j)+10^-25);
un(i,j)=(anv(i,j)-p1)*lpr(i,j);
ue(i,j)=(aev(i,j)-p2)*lpr(i,j);
end
end
for l=1:3
resl=0;
for i=2:ny
for j=2:nx+1
res(i,j)=sv(i,j)-(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+apv(i,j)*v(i,j,1));
resl=resl+abs(res(i,j));
res(i,j)=(res(i,j)-ls(i,j)*res(i-1,j)-lw(i,j)*res(i,j-1))*lpr(i,j);
end
end
if l==1
resor=resl;
end
% rsm=resl/(resor+10^-25)
resl;
for i=(ny):-1:2
for j=(nx+1):-1:2
res(i,j)=res(i,j)-un(i,j)*res(i+1,j)-ue(i,j)*res(i,j+1);
v(i,j,1)=v(i,j,1)+res(i,j);
end
end
end
for i=1:ny
for j=1:nx
%mp(i,j)=raw*dx*((u(i+1,j+1,1)-u(i+1,j,1))+(v(i+1,j+1,1)-v(i,j+1,1)));
end
end
%norm(mp);%mp become nonzero!
end
%%%%%%%%%%%%%%%%%%%pressure correction%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%co. for pressure correction
for i=2:ny+1
for j=2:nx
if i==2
aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
asu(i,j)=0;%d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
%awu(:,2)=0;
%aeu(:,nx)=0;
apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j)+0.6*d x/u(i,j,1);%(dy*(1.2*nx/ny)/nx+0.1*(u(i+1,j,1)-u(i,j,1)))/u(i,j,1);
else if i==ny+1
aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
anu(i,j)=0;%d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
%awu(:,2)=0;
%aeu(:,nx)=0;
apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j)+0.6*d x/u(i,j,1);%+(dy*(1.2*nx/ny)/nx+0.1*(u(i-1,j,1)-u(i,j,1)))/u(i,j,1);
else
aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
%awu(:,2)=0;
%aeu(:,nx)=0;
apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j);
end
end
%u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i-1,j-1,1)-p(i-1,j,1))*dy)/apu(i,j);
end
end
%u(:,nx,1)=u(:,nx,1)*(sum(u(:,1,1))/sum(u(:,nx,1)));
%u(:,nx+1,1)=u(:,nx,1);
for i=2:ny
for j=2:nx+1
if j==2
aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
awv(i,j)=0;%2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j-1,1)+0.25*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
%v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+(p(i-1,j-1,1)-p(i,j-1,1))*dy)/apv(i,j);
else if j==nx+1
aev(i,j)=0;%2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j,1)+0.25*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
%v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+(p(i-1,j-1,1)-p(i,j-1,1))*dy)/apv(i,j);
else
aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
%v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+(p(i-1,j-1,1)-p(i,j-1,1))*dy)/apv(i,j);
end
end
end
end
for i=1:ny
for j=2:nx-1
aepc(i,j)=raw*dy^2/apu(i+1,j+1);
awpc(i,j)=raw*dy^2/apu(i+1,j);
end
end
for i=2:ny-1
for j=1:nx
anpc(i,j)=raw*dx^2/apv(i+1,j+1);
aspc(i,j)=raw*dx^2/apv(i,j+1);
end
end
for i=1:ny
awpc(i,1)=0;
aepc(i,1)=raw*dy^2/apu(i+1,2);
awpc(i,nx)=raw*dy^2/apu(i+1,nx);
aepc(i,nx)=0;
end
for j=1:nx
anpc(1,j)=raw*dx^2/apv(2,j+1);
aspc(1,j)=0;
anpc(ny,j)=0;
aspc(ny,j)=raw*dx^2/apv(ny,j+1);
end
aepc(15:16,39)=0;
awpc(15:16,42)=0;
anpc(14,40:41)=0;
aspc(17,40:41)=0;
for i=1:ny
for j=1:nx
appc(i,j)=aepc(i,j)+awpc(i,j)+anpc(i,j)+aspc(i,j);
end
end
for i=1:ny
for j=1:nx
mp(i,j)=raw*dx*((u(i+1,j+1,1)-u(i+1,j,1))+(v(i+1,j+1,1)-v(i,j+1,1)));
end
end


%SIP METHOD for pressure correction
sp=zeros(ny+2,nx+2);
urf=0.02;
%co. for SIP method
ls=zeros(ny+2,nx+2);
lw=zeros(ny+2,nx+2);
un=zeros(ny+2,nx+2);
ue=zeros(ny+2,nx+2);
lpr=zeros(ny+2,nx+2);
res=zeros(ny+2,nx+2);
resl=0;
for i=2:ny+1
for j=2:nx+1
appc(i-1,j-1)=-appc(i-1,j-1)/0.2;
sp(i,j)=0.8*appc(i-1,j-1)*pc(i,j,1)+mp(i-1,j-1);
end
end
for j=2:nx+1
for i=2:ny+1
lw(i,j)=awpc(i-1,j-1)/(1+0.9*un(i,j-1));
ls(i,j)=aspc(i-1,j-1)/(1+0.9*ue(i-1,j));
p1=0.9*lw(i,j)*un(i,j-1);
p2=0.9*ls(i,j)*ue(i-1,j);
lpr(i,j)=1/(appc(i-1,j-1)+p1+p2-lw(i,j)*ue(i,j-1)-ls(i,j)*un(i-1,j)+10^-25);
un(i,j)=(anpc(i-1,j-1)-p1)*lpr(i,j);
ue(i,j)=(aepc(i-1,j-1)-p2)*lpr(i,j);
end
end
for l=1:5
resl=0;
for i=2:ny+1
for j=2:nx+1
res(i,j)=sp(i,j)-(aepc(i-1,j-1)*pc(i,j+1,1)+awpc(i-1,j-1)*pc(i,j-1,1)+anpc(i-1,j-1)*pc(i+1,j,1)+aspc(i-1,j-1)*pc(i-1,j,1)+appc(i-1,j-1)*pc(i,j,1));
resl=resl+abs(res(i,j));
res(i,j)=(res(i,j)-ls(i,j)*res(i-1,j)-lw(i,j)*res(i,j-1))*lpr(i,j);
end
end
if l==1
resor=resl;
end
% rsm=resl/(resor+10^-25)
resl;
for i=(ny+1):-1:2
for j=(nx+1):-1:2
res(i,j)=res(i,j)-un(i,j)*res(i+1,j)-ue(i,j)*res(i,j+1);
pc(i,j,1)=pc(i,j,1)+res(i,j);
end
end
end

%keep the corner pressure tobe 0
ref=pc(2,2,1);
for i=2:ny+1
for j=2:nx+1
pc(i,j,1)=pc(i,j,1)-ref;
end
end
%correct the pressure field
for ll=1:200
for i=2:ny+1
for j=2:nx+1
% pc(i,j,1)=(pc(i,j+1,1)*aepc(i-1,j-1)+pc(i,j-1,1)*awpc(i-1,j-1)+pc(i+1,j,1)*anpc(i-1,j-1)+pc(i-1,j,1)*aspc(i-1,j-1)-mp(i-1,j-1))/appc(i-1,j-1);
end
end
end

for i=1:ny
for j=1:nx
p(i,j,1)=p(i,j,1)+alpha*pc(i+1,j+1,1);
end
end
%correct velocity field
for i=2:ny+1
for j=2:nx
u(i,j,1)=u(i,j,1)+dx/apu(i,j)*(pc(i,j,1)-pc(i,j+1,1));
end
end
for i=2:ny
for j=2:nx
v(i,j,1)=v(i,j,1)+dy/apv(i,j)*(pc(i,j,1)-pc(i+1,j,1));
end
end
%check for residual
err=norm(mp(:,)
errdata(q)=err;
u(16:17,40:42)=0;
v(15:17,41:42)=0;
if err<10^-8
break
end
clf
subplot(2,2,1),contourf(u(:,:,1),30);axis equal
subplot(2,2,2), plot(errdata(1,1:q));
subplot(2,2,3),contourf(pc(:,:,1),30);axis equal
subplot(2,2,4),contourf(p(:,:,1),30);axis equal
hold off;
drawnow
end

[prasanthnitt]

Please specify the problem that you have solved and the algorithm that you have used for the same.

[mgedwards]

Best explain the problem
You have left some bugs

[kalanpi4]

Hey,

Can you explain how have you set up the grid? I'm having problems implementing SIMPLE in MATLAB for flow between plates.

Thanks in advance,

Kalyan