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

HELP!Problems on SIMPLE method

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

Reply
 
LinkBack Thread Tools Display Modes
Old   July 8, 2011, 07:15
Default HELP!Problems on SIMPLE method
  #1
Member
 
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 6
houkensjtu is on a distinguished road
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
houkensjtu is offline   Reply With Quote

Old   July 8, 2011, 11:22
Default
  #2
Member
 
SM
Join Date: Dec 2010
Posts: 61
Rep Power: 6
canopus is on a distinguished road
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
canopus is offline   Reply With Quote

Old   July 9, 2011, 10:49
Default
  #3
Member
 
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 6
houkensjtu is on a distinguished road
Quote:
Originally Posted by canopus View Post
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?
houkensjtu is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
SIMPLE family and Fractional Step Method Geon-Hong Main CFD Forum 2 April 5, 2010 18:05
SIMPLE method for inviscid flow abcdef123 Main CFD Forum 0 February 26, 2010 09:24
Apply Multigrid method to SIMPLE Algorithm yang Main CFD Forum 1 February 25, 2006 12:28
need help in fotran code for simple method saritha Main CFD Forum 3 April 30, 2004 12:55
SIMPLE method for compressible flow in pipe Reza Main CFD Forum 0 August 27, 2002 22:43


All times are GMT -4. The time now is 04:57.