|
[Sponsors] |
July 8, 2011, 07:15 |
HELP!Problems on SIMPLE method
|
#1 |
Member
HouKen
Join Date: Jul 2011
Posts: 67
Rep Power: 14 |
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: 97
Rep Power: 15 |
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: 67
Rep Power: 14 |
Quote:
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 | Search this Thread |
Display Modes | |
|
|
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 08:24 |
Apply Multigrid method to SIMPLE Algorithm | yang | Main CFD Forum | 1 | February 25, 2006 11: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 |