# HELP!Problems on SIMPLE method

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

 July 8, 2011, 07:15 HELP!Problems on SIMPLE method #1 Member   HouKen Join Date: Jul 2011 Posts: 66 Rep Power: 7 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

 July 8, 2011, 11:22 #2 Member   SM Join Date: Dec 2010 Posts: 90 Rep Power: 8 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

July 9, 2011, 10:49
#3
Member

HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 7
Quote:
 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?

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Geon-Hong Main CFD Forum 2 April 5, 2010 18:05 abcdef123 Main CFD Forum 0 February 26, 2010 09:24 yang Main CFD Forum 1 February 25, 2006 12:28 saritha Main CFD Forum 3 April 30, 2004 12:55 Reza Main CFD Forum 0 August 27, 2002 22:43

All times are GMT -4. The time now is 23:19.