[houkensjtu]

I'm a master student and trying to do a laminar simulation by SIMPLE method.
I use matlab and the result seems not good.
The most serious problem I think is,It seems very difficult to get convergence solution on Pressure correction.(I have tried Jacobi and ADI)NOT ONLY the absolute value of pressure grows with iteration，but ALSO the pressure drop grows,while the contour shape of pressure drop seems NOT strange.

The code is writing by matlab,based on staggered grid by Patankar.(I didn't do residual check in the program just for easy to read.The fact is even I check for residual,it can't get a convergence result)

clear all;
%initialization
nx=30;
ny=10;
%velcity field
u=zeros(ny,nx+1,3);
v=zeros(ny+1,nx,3);
%pressure
p=zeros(ny,nx,3);
%co. for momentum equation
aeu=zeros(ny,nx);
awu=zeros(ny,nx);
asu=zeros(ny,nx);
anu=zeros(ny,nx);
apu=zeros(ny,nx);
aev=zeros(ny,nx);
awv=zeros(ny,nx);
asv=zeros(ny,nx);
anv=zeros(ny,nx);
apv=zeros(ny,nx);
%co. for pressure correction equation
pce=zeros(ny,nx);
pcw=zeros(ny,nx);
pcs=zeros(ny,nx);
pcn=zeros(ny,nx);
aep=zeros(ny,nx);
awp=zeros(ny,nx);
asp=zeros(ny,nx);
anp=zeros(ny,nx);
app=zeros(ny,nx);
pc=zeros(ny,nx,3);
%mass unbalance
mp=zeros(ny,nx);
errp=zeros(ny,nx);
raw=1;
dx=0.1;
dy=0.1;
d=1;
u(2:ny-1,:,=30;
for q=1:5000
for ll=1:30
for i=2:ny-1
for j=2:nx
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,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)));
anu(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)));
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,j-1,1)-p(i,j,1))*dy)/apu(i,j);
end
end
for i=2:ny
for j=2:nx-1
aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i-1,j+1,1)+0.5*u(i,j+1,1))/d))^5)+max(0,-raw*dy*(0.5*u(i-1,j+1,1)+0.5*u(i,j+1,1)));
awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i-1,j,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i-1,j,1)+0.5*u(i,j,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)-p(i,j,1))*dy)/apv(i,j);
end
end
end


for ll=1:100
for i=2:ny-1
for j=2:nx-1
cep(i,j)=raw*dy*(0.375*(u(i,j,1)+u(i,j+1,1))+(u(i-1,j,1)+u(i+1,j,1)+u(i-1,j+1,1)+u(i+1,j+1,1))/16);
cwp(i,j)=raw*dy*(0.375*(u(i,j,1)+u(i,j-1,1))+(u(i-1,j,1)+u(i+1,j,1)+u(i-1,j-1,1)+u(i+1,j-1,1))/16);
cnp(i,j)=raw*dx*(0.375*(v(i+1,j,1)+v(i,j,1))+(v(i, j-1,1)+v(i,j+1,1)+v(i+1,j-1,1)+v(i+1,j+1,1))/16);
csp(i,j)=raw*dx*(0.375*(v(i-1,j,1)+v(i,j,1))+(v(i,j-1,1)+v(i,j+1,1)+v(i-1,j-1,1)+v(i-1,j+1,1))/16);
pcep(i,j)=cep(i,j)/d;
pcwp(i,j)=cwp(i,j)/d;
pcnp(i,j)=cnp(i,j)/d;
pcsp(i,j)=csp(i,j)/d;
aep(i,j)=d*(max(0,(1-0.1*abs(pcep(i,j)))^5)+max(-pcep(i,j),0));
awp(i,j)=d*(max(0,(1-0.1*abs(pcwp(i,j)))^5)+max(pcwp(i,j),0));
anp(i,j)=d*(max(0,(1-0.1*abs(pcnp(i,j)))^5)+max(-pcnp(i,j),0));
asp(i,j)=d*(max(0,(1-0.1*abs(pcsp(i,j)))^5)+max(pcsp(i,j),0));
app(i,j)=aep(i,j)+awp(i,j)+asp(i,j)+anp(i,j);
end
end
for i=2:ny-1
for j=2:nx-1
aepc(i,j)=2*raw*dy^2/(app(i,j)+app(i,j+1));
awpc(i,j)=2*raw*dy^2/(app(i,j)+app(i,j-1));
anpc(i,j)=2*raw*dy^2/(app(i,j)+app(i+1,j));
aspc(i,j)=2*raw*dy^2/(app(i,j)+app(i-1,j));
end
end
for i=2:ny-1
for j=2:nx-1
mp(i,j)=raw*dx*((u(i,j+1,1)-u(i,j-1,1))/2+(v(i+1,j,1)-v(i-1,j,1))/2);
end
end

for i=2:ny-1
aepc(i,nx-1)=raw*dy^2/app(i,nx-1);
end
for i=2:ny-1
awpc(i,2)=raw*dy^2/app(i,2);
end
for j=2:nx-1
anpc(ny-1,j)=raw*dy^2/app(ny-1,j);
end
for j=2:nx-1
aspc(2,j)=raw*dy^2/app(2,j);
end

for i=2:ny-1
for j=2:nx-1
appc(i,j)=aepc(i,j)+awpc(i,j)+anpc(i,j)+aspc(i,j);
end
end

for i=2:ny-1
for j=2:nx-1
pc(i,j,1)=(pc(i,j+1,1)*aepc(i,j)+pc(i,j-1,1)*awpc(i,j)+pc(i+1,j,1)*anpc(i,j)+pc(i-1,j,1)*aspc(i,j)-mp(i,j))/appc(i,j);
end
end
end





for i=2:ny-1
pc(i,1,1)=pc(i,2,1);
pc(i,nx,1)=pc(i,nx-1,1);
end
for j=2:nx-1
pc(1,j,1)=pc(2,j,1);
pc(ny,j,1)=pc(ny-1,j,1);
end
pc(1,1,1)=(pc(1,2,1)+pc(2,1,1))/2;
pc(1,nx,1)=(pc(1,nx-1,1)+pc(2,nx,1))/2;
pc(ny,1,1)=0.5*(pc(ny,2,1)+pc(ny-1,1,1));
pc(ny,nx,1)=0.5*(pc(ny-1,nx,1)+pc(ny,nx-1,1));

p(:,:,1)=p(:,:,1)+0.01*pc(:,:,1);
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,1)-pc(i,j,1));
end
end
for i=2:ny
for j=2:nx-1
v(i,j,1)=v(i,j,1)+dy/apv(i,j)*(pc(i-1,j,1)-pc(i,j,1));
end
end
err=norm(mp(:,)
if err<10^-10
break
end

end

[canopus]

Hello!
If you post a long piece of code and expect someone to debug it then you should be really lucky!

maybe you can look into this -

http://cfd.ce.gatech.edu/docs/CEE7751_7.pdf

http://www-math.mit.edu/cse/codes/mi...vierstokes.pdf

[houkensjtu]

Originally Posted by canopus

Hello!
If you post a long piece of code and expect someone to debug it then you should be really lucky!

maybe you can look into this -

http://cfd.ce.gatech.edu/docs/CEE7751_7.pdf

http://www-math.mit.edu/cse/codes/mi...vierstokes.pdf

Thank u very much!
I post another question when I read codes.If u are familiar with SIMPLE would u please also answer that problem?