
[Sponsors] 
August 20, 2012, 03:10 

#21 
Senior Member
Join Date: Aug 2011
Posts: 271
Rep Power: 9 

August 20, 2012, 03:37 

#22  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 723
Rep Power: 19 
Quote:


August 20, 2012, 04:04 

#23  
Senior Member
Join Date: Aug 2011
Posts: 271
Rep Power: 9 
Quote:
His linear system is fine because it works in artificial compressibility case for the velocity componnents. Recall indeed that in artificial compressibility you don't have a linear system to solve for pressure since it is determined explicitly with an expression like P(i,j)_k+1= P(i,j)_k + c*div(U*). k is the indice loop to reach the steady state, c is a pseudo celerity of pressure waves. So if his SIMPLE does not work it should no be due to the linear system solver but rather due to its implementation which should have a bugg. 

August 20, 2012, 04:08 

#24 
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 723
Rep Power: 19 

August 20, 2012, 17:55 

#25 
Member
Junaid Ahmad Khan
Join Date: Mar 2010
Location: Islamabad
Posts: 33
Rep Power: 9 
Hi Every one,
Thanks for all your support tonight i have seen results that i was hopping to see for about more then a week. and some how i have solve the Lid Driven Cavity by simple algorithm. I will post the code and procedure later for future reference. well there is still order of convergence problems and the understanding of relaxation factor which will be discuss later. Thanks Again. Regards Junaid 

August 20, 2012, 19:27 

#26 
Member
Michael Moor
Join Date: May 2012
Location: Ireland
Posts: 30
Rep Power: 7 
Hi Junaid,
Congratulations on your results!! Please may you also send the program on to me? I am curious as to how you implemented the boundary conditions for the vertical walls. 

August 21, 2012, 01:55 

#27  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 723
Rep Power: 19 
Quote:
In case of SIMPLE algo enclosed with walls (lid driven cavity) your pressure equation shall be all neumann BC type problem, which should be singular. How did you solve that. Probably that is main issue for your convergence problems. 

August 21, 2012, 03:25 

#28 
Member
Ren/Xingyue
Join Date: Jan 2010
Location: Nagoya , Japan
Posts: 44
Rep Power: 9 
good topic!
I think artificial compressiblity method is more stable then SIMPLE with a bad linear solver in your case. Junaid, can you tell me your solver? Maybe you didn't use a suitable solver. Last edited by hilllike; August 21, 2012 at 03:59. 

August 21, 2012, 04:41 

#29 
Senior Member
Join Date: Aug 2011
Posts: 271
Rep Power: 9 
As I mentioned it previously there is no linear system to solve for pressure in the artificial compressibility method.


August 21, 2012, 04:55 

#30  
Member
Ren/Xingyue
Join Date: Jan 2010
Location: Nagoya , Japan
Posts: 44
Rep Power: 9 
Quote:
I don't know if his SIMPLE code is two phase flow model. If it is the linear solver may be the problem. I talked about the solver in SIMPLE code. 

August 21, 2012, 05:54 

#31 
Member
Junaid Ahmad Khan
Join Date: Mar 2010
Location: Islamabad
Posts: 33
Rep Power: 9 
The problem with convergence is that when i let the program to iterate to go to more then about 8,000 the streamlines become distorted. As for pressure i solve it explicitly. as given below
Kindly go through it and try to debug it. it is a working code #include <iostream> #include <vector> #include <fstream> #include <math.h> #include <time.h> #include <algorithm> #include <conio.h> #include <string> #include <stdlib.h> using namespace std; //Number of NODEs #define NODE 81 //Write file After that number of Iterations #define AutoSave 1000 //UnderRelaxation Factor for p,u,v float urfp=0.0001; float urfu=0.3; float urfv=0.7; // //Boundary conditions velocitys float B[2]={1,0}; // //X and Y coordinate float x[NODE] = {0}; float y[NODE] = {0}; // //Uvelocity, Vvelocity, and Pressure float u[NODE][NODE+1] = {0}; float us[NODE][NODE+1] = {0}; float du[NODE][NODE+1] = {0}; float uo[NODE][NODE+1] = {0}; float v[NODE+1][NODE] = {0}; float vs[NODE+1][NODE] = {0}; float dv[NODE+1][NODE] = {0}; float vo[NODE+1][NODE] = {0}; float p[NODE+1][NODE+1] = {0.0}; float ps[NODE+1][NODE+1] = {0.0}; float pp[NODE+1][NODE+1] = {0.0}; float po[NODE+1][NODE+1] = {0.0}; // void TDMsolve (int n, float *a, float *b, float *c, float *v, float *x) { /** * n  number of equations * a  subdiagonal (means it is the diagonal below the main diagonal)  indexed from 1..n1 * b  the main diagonal * c  supdiagonal (means it is the diagonal above the main diagonal)  indexed from 0..n2 * v  right part * x  the answer */ for (int i = 1; i < n; i++) { double m = a[i]/b[i1]; b[i] = b[i]  m*c[i1]; v[i] = v[i]  m*v[i1]; } x[n1] = v[n1]/b[n1]; for (int i = n  2; i >= 0; i) x[i]=(v[i]c[i]*x[i+1])/b[i]; } float MAX(float f,float d) { if(f>d) return f; if(d>f) return d; else return f; } void writetoFile(int Nx,int Ny,float comp_time,float er, char* name) { ofstream outfile2; outfile2.open(name); outfile2<<"TITLE = \"Computation Time = "<<comp_time<<" ms, Error = "<<er<<"\"\nVARIABLES = \"X\", \"Y\", \"u(x,y,t)\", \"v(x,y,t)\", \"p(x,y,t)\"\n\n"; outfile2<<"ZONE I="<<Nx<<", J="<<Ny<<", ZONETYPE=Ordered, DATAPACKING=POINT\n"; for (int row=0; row<Ny; row++) { for (int col=0; col<Nx; col++) { outfile2<<x[col]<<"\t"<<y[row]<<"\t" <<(0.5*(u[col][row]+u[col][row+1]))<<"\t" <<(0.5*(v[col][row]+v[col+1][row]))<<"\t" <<(0.25*(p[col][row]+p[col+1][row]+p[col][row+1]+p[col+1][row+1]))<<"\n"; } } outfile2.close(); } //Ser Boundary Conditions void setBC(int Nx,int Ny) { for(int row=0;row<Ny;row++) { v[0][row] = 2.0*B[1]v[1][row];//West Wall v[Ny][row] = 2.0*B[1]v[Nx1][row];//East wall } for(int col=0;col<Nx;col++) { u[col][0]=2.0*B[0]u[col][1];//North wall u[col][Nx]=2.0*B[1]u[col][Nx1];//South wall } } void copyValues(int Nx,int Ny) { for (int col=0; col<Nx; col++) { for (int row=0; row<Ny; row++) { ps[col][row] = p[col][row]; us[col][row] = u[col][row]; vs[col][row] = v[col][row]; } us[col][Ny] = u[col][Ny]; vs[Nx][col] = v[Nx][col]; ps[Nx][col] = p[Nx][col]; ps[col][Ny] = p[col][Ny]; } ps[Nx][Ny] = p[Nx][Ny]; } void CPO(int Nx,int Ny) { //Copy Old Values for (int col=0; col<Nx; col++) { for (int row=0; row<Ny; row++) { po[col][row] = p[col][row]; //uo[col][row] = us[col][row]; //vo[col][row] = vs[col][row]; } //uo[col][Ny] = us[col][Ny]; //vo[Nx][col] = vs[Nx][col]; po[Nx][col] = p[Nx][col]; po[col][Ny] = p[col][Ny]; } po[Nx][Ny] = p[Nx][Ny]; } /***********Calculate Pressure Correction *************/ void calcPp(int Nx,int Ny,float A,float G) { float aW=0.0,aE=0.0,aS=0.0,aN=0.0,aP=0.0; float bpw=0.0,bpe=0.0,bps=0.0,bpn=0.0,bp=0.0; for(int row=1; row<Ny; row++) { for(int col=1; col<Nx; col++) { aW=du[col1][row]*A; aE=du[col][row]*A; aS=dv[col][row]*A; aN=dv[col][row1]*A; aP=aW+aE+aS+aN; bpw=us[col1][row]*A; bpe=us[col][row]*A; bps=vs[col][row]*A; bpn=vs[col][row1]*A; bp=bpwbpe+bpsbpn; pp[col][row]=( aW*pp[col1][row] + aE*pp[col+1][row] + aS*pp[col][row+1] + aN*pp[col][row1] + bp )/aP; //a(I,J) p'(I,J)=a(I+1,J) x p'(I+1,J ) + a(I1,J) x p'(I1,J ) + a(I,J+1) x p'(I,J+1) + a(I,J1) x p'(I,J1) + b'(I,J) } } } void PressureCorrection(int Nx,int Ny) { for (int col=1; col<Nx; col++) { for (int row=1; row<Ny; row++) { p[col][row] = ps[col][row]+(urfp*pp[col][row]); // p^new = p* + alpha_p x p' } } } void UCorrection(int Nx,int Ny,float A,float G) { //UVelocity Corrector // float d=0.0; for(int row=1; row<Ny; row++) { for(int col=1; col<Nx1; col++) { d=du[col][row]; u[col][row]=us[col][row]+(d*(pp[col][row]pp[col+1][row])); //I have used farword Stagered in U // u = u* + d(I,J) (p'(I,J)p'(I+1,J)) } } } //Uunder relaxation void U_urf(int Nx,int Ny,float A,float G) { //UVelocity Corrector float Fw=0.0,Fe=0.0,Fs=0.0,Fn=0.0,dF=0.0; float FW=0.0,FE=0.0,FP=0.0,FSE=0.0,FSW=0.0,FNE=0.0,FNW=0 .0; float D=0.0; float dxw=0.0,dxe=0.0,dyn=0.0,dys=0.0; float aw=0.0,ae=0.0,as=0.0,an=0.0,ap=0.0,apo=0.0; float U=0.0; float a1[NODE]={0},b1[NODE]={0},c1[NODE]={0},v1[NODE]={0},x1[NODE]={0}; int sb=0; for(int row=1; row<Ny; row++) { for(int col=0; col<Nx; col++) { if(col==0col==(Nx1)) { a1[sb]=0.0; c1[sb]=0.0; b1[sb]=1.0; v1[sb]=u[col][row]; } else if(col>0&&col<(Nx1)) { FP=u[col][row]; FW=u[col1][row]; FE=u[col+1][row]; FSE=v[col+1][row]; FSW=v[col][row]; FNE=v[col+1][row1]; FNW=v[col][row1]; Fw=A*(FP+FW)/2.0; Fe=A*(FE+FP)/2.0; Fs=A*(FSE+FSW)/2.0; Fn=A*(FNE+FNW)/2.0; dF=FeFw+FnFs; aw=MAX(G+Fw/2.0,MAX(Fw,0.0)); ae=MAX(GFe/2.0,MAX(Fe,0.0)); as=MAX(G+Fs/2.0,MAX(Fs,0.0)); an=MAX(GFn/2.0,MAX(Fn,0.0)); ap=(aw+ae+as+an+dF)/urfu; a1[sb]=aw; c1[sb]=ae; b1[sb]=ap; v1[sb]=( as*u[col][row+1] + an*u[col][row1] + A*(p[col][row]p[col+1][row]) + ((1urfu)*ap)*uo[col][row] ); //a(i,j)/alpha_u x u(i,j) = Sigma (a_nb x v_nb) + (p(I,J)+p(I+1,J)) x A +[(1alpha_u)a(i,j)/alpha_u] x u^(n1)_(I,j) } sb++; } TDMsolve (NODE, a1, b1, c1, v1, x1 ); for(int k=0;k<Ny;k++){ u[k][row]=x1[k]; } sb=0; } } void VCorrection(int Nx,int Ny,float A,float G) { //VVelocity Corrector float v_calc=0.0,d=0.0; for(int col=1; col<Nx; col++) { for(int row=1; row<Ny1; row++) { d=dv[col][row]; v[col][row]=vs[col][row]+(d*(pp[col][row+1]pp[col][row])); //From top to bottom J increase // v = v* + d(I,J) (p'(I,J+1)p'(I,J)) } } } //Vunder relaxation void V_urf(int Nx,int Ny,float A,float G) { float Fw=0.0,Fe=0.0,Fs=0.0,Fn=0.0,dF=0.0; float FN=0.0,FS=0.0,FP=0.0,FSE=0.0,FSW=0.0,FNE=0.0,FNW=0 .0; float D=0.0; float dxw=0.0,dxe=0.0,dyn=0.0,dys=0.0; float aw=0.0,ae=0.0,as=0.0,an=0.0,ap=0.0,apo=0.0; float a2[NODE]={0},b2[NODE]={0},c2[NODE]={0},v2[NODE]={0},x2[NODE]={0}; int sb=0; for(int col=1; col<Nx; col++) { for(int row=0; row<Ny; row++) { if(row==0row==(Ny1)) { a2[sb]=0.0; c2[sb]=0.0; b2[sb]=1.0; v2[sb]=v[col][row]; } else if(row>0&&row<(Ny1)) { FSW=(u[col1][row+1]); FNW=(u[col1][row]); FNE=(u[col][row]); FSE=(u[col][row+1]); FP=(v[col][row]); FS=(v[col][row+1]); FN=(v[col][row1]); Fw=A*(FSW+FNW)/2.0; Fe=A*(FNE+FSE)/2.0; Fs=A*(FS+FP)/2.0; Fn=A*(FN+FP)/2.0; dF=FeFw+FnFs; D=G; aw=MAX(D+Fw/2.0,MAX(Fw,0.0)); ae=MAX(DFe/2.0,MAX(Fe,0.0)); as=MAX(D+Fs/2.0,MAX(Fs,0.0)); an=MAX(DFn/2.0,MAX(Fn,0.0)); ap=(aw+ae+as+an+dF)/urfv; a2[sb]=an; c2[sb]=as; b2[sb]=ap; v2[sb]=( aw*v[col1][row] + ae*v[col+1][row] + A*(p[col][row+1]p[col][row]) + ((1urfv)*ap)*vo[col][row] //a(I,j)/alpha_v x v(I,j) = Sigma (a_nb x v_nb) + (p(I,J+1)+p(I,J)) x A +[(1alpha_v)a(I,j)/alpha_v] x v^(n1)_(I,j) ); } sb++; } TDMsolve (NODE, a2, b2, c2, v2, x2); for(int k=0;k<Ny;k++){ v[col][k]=x2[k]; } sb=0; } } void setWallGradient(int Nx,int Ny) { /***********setWallPressGradientToZero************* *************/ for (int col=0; col<=Nx; col++) { //North p[col][0] = p[col][1]; //South p[col][Ny] = p[col][Ny1]; } for (int row=0; row<=Ny; row++) { //West p[0][row] = p[1][row]; //East p[Nx][row] = p[Nx1][row]; } } void calcUs(int Nx,int Ny,float A,float G) { /***********Calculate uVelocity at Wall *************/ float Fw=0.0,Fe=0.0,Fs=0.0,Fn=0.0,dF=0.0; float FW=0.0,FE=0.0,FP=0.0,FSE=0.0,FSW=0.0,FNE=0.0,FNW=0 .0; float D=0.0; float dxw=0.0,dxe=0.0,dyn=0.0,dys=0.0; float aw=0.0,ae=0.0,as=0.0,an=0.0,ap=0.0,apo=0.0; float U=0.0; float a1[NODE]={0},b1[NODE]={0},c1[NODE]={0},v1[NODE]={0},x1[NODE]={0}; int sb=0; for(int row=1; row<Ny; row++) { for(int col=0; col<Nx; col++) { if(col==0col==(Nx1)) { a1[sb]=0.0; c1[sb]=0.0; b1[sb]=1.0; v1[sb]=us[col][row]; du[col][row]=A; } else if(col>0&&col<(Nx1)) { FP=us[col][row]; FW=us[col1][row]; FE=us[col+1][row]; FSE=vs[col+1][row]; FSW=vs[col][row]; FNE=vs[col+1][row1]; FNW=vs[col][row1]; Fw=A*(FP+FW)/2.0; Fe=A*(FE+FP)/2.0; Fs=A*(FSE+FSW)/2.0; Fn=A*(FNE+FNW)/2.0; dF=FeFw+FnFs; aw=MAX(G+Fw/2.0,MAX(Fw,0.0)); ae=MAX(GFe/2.0,MAX(Fe,0.0)); as=MAX(G+Fs/2.0,MAX(Fs,0.0)); an=MAX(GFn/2.0,MAX(Fn,0.0)); ap=aw+ae+as+an+dF; du[col][row]=A/ap; a1[sb]=aw; c1[sb]=ae; b1[sb]=ap; v1[sb]=( as*us[col][row+1] + an*us[col][row1] + A*(ps[col][row]ps[col+1][row]) ); //a(I,j) x u*(I,j) = Sigma (a_nb x u*_nb) + (p*(I,J)+p*(I+1,J)) x A } sb++; } TDMsolve (NODE, a1, b1, c1, v1, x1 ); for(int k=0;k<Ny;k++){ us[k][row]=x1[k]; } sb=0; } } void calcVs(int Nx,int Ny,float A,float G) { /***********Calculate vVelocity at Wall *************/ float Fw=0.0,Fe=0.0,Fs=0.0,Fn=0.0,dF=0.0; float FN=0.0,FS=0.0,FP=0.0,FSE=0.0,FSW=0.0,FNE=0.0,FNW=0 .0; float D=0.0; float dxw=0.0,dxe=0.0,dyn=0.0,dys=0.0; float aw=0.0,ae=0.0,as=0.0,an=0.0,ap=0.0,apo=0.0; float a2[NODE]={0},b2[NODE]={0},c2[NODE]={0},v2[NODE]={0},x2[NODE]={0}; int sb=0; for(int col=1; col<Nx; col++) { for(int row=0; row<Ny; row++) { if(row==0row==(Ny1)) { a2[sb]=0.0; c2[sb]=0.0; b2[sb]=1.0; v2[sb]=vs[col][row]; dv[col][row]=A; } else if(row>0&&row<(Ny1)) { FSW=(us[col1][row+1]); FNW=(us[col1][row]); FNE=(us[col][row]); FSE=(us[col][row+1]); FP=(vs[col][row]); FS=(vs[col][row+1]); FN=(vs[col][row1]); Fw=A*(FSW+FNW)/2.0; Fe=A*(FNE+FSE)/2.0; Fs=A*(FS+FP)/2.0; Fn=A*(FN+FP)/2.0; dF=FeFw+FnFs; D=G; aw=MAX(D+Fw/2.0,MAX(Fw,0.0)); ae=MAX(DFe/2.0,MAX(Fe,0.0)); as=MAX(D+Fs/2.0,MAX(Fs,0.0)); an=MAX(DFn/2.0,MAX(Fn,0.0)); ap=aw+ae+as+an+dF; dv[col][row]=A/ap; a2[sb]=an; c2[sb]=as; b2[sb]=ap; v2[sb]=( aw*vs[col1][row] + ae*vs[col+1][row] + A*(ps[col][row+1]ps[col][row]) ); //a(I,j) x v*(I,j) = Sigma (a_nb x v*_nb) + (p*(I,J+1)+p*(I,J)) x A } sb++; } TDMsolve (NODE, a2, b2, c2, v2, x2); for(int k=0;k<Ny;k++){ vs[col][k]=x2[k]; } sb=0; } } float calcError(int Nx,int Ny) { float temp = fabs(p[0][0]  po[0][0]); float val; for (int row=0; row<Ny+1; row++) for (int col=0; col<Nx+1; col++) { val = fabs(p[col][row]  po[col][row]); if (val>temp) temp = val; } return temp; } /***********MAIN*************/ int main() { /**********Variables and Constants******/ float Lx = 1.0; float Ly = 1.0; int nx = NODE; int ny = NODE; float deltaX = Lx/(nx1); float deltaY = Ly/(ny1); float RE = 1000.0; float Gm = (1.0/RE); float er = 100.0; float tempx = 0.0; for (int col=0; col<nx; col++) {x[col] = tempx; tempx +=deltaX;} float tempy = 1.0; for (int row=0; row<ny; row++) {y[row] = tempy; tempy =deltaY;} float Area=deltaX; //***********Display Grid*************/ //Initialize Grid for (int col=0; col<nx; col++) { for (int row=0; row<ny; row++) { p[col][row] = 0.1; pp[col][row] = 0.0; u[col][row] = 0.0; v[col][row] = 0.0; } u[col][ny] = 0.0; v[nx][col] = 0.0; p[nx][col] = 0.1; p[col][ny] = 0.1; pp[nx][col] = 0.0; pp[col][ny] = 0.0; } p[nx][ny] = 0.1; pp[nx][ny] = 0.0; // setBC(nx,ny); copyValues(nx,ny); CPO(nx,ny); // time_t end, start_t; double time; start_t = clock(); int count = 1; int it=0; float it1=0; char nm[20]; float Er=0.000001; //File Name itoa(it,nm,10); int ch=0; for(;nm[ch]!='\0';ch++){} nm[ch]='.'; nm[ch+1]='d'; nm[ch+2]='a'; nm[ch+3]='t'; nm[ch+4]='\0'; writetoFile(nx,ny,0,0,nm); cout<<"\nAutosave"<<endl; int iter=0; // //Iteration Loop while (er>=Er) { //Set Boundary Condition calcUs(nx,ny,Area,Gm); calcVs(nx,ny,Area,Gm); //Calculate p' calcPp(nx,ny,Area,Gm); //Corrected Pressure PressureCorrection(nx,ny); //Set Wall Gradient setWallGradient(nx,ny); UCorrection(nx,ny,Area,Gm); VCorrection(nx,ny,Area,Gm); setBC(nx,ny); U_urf(nx,ny,Area,Gm); V_urf(nx,ny,Area,Gm); setBC(nx,ny); //Copy Values From Field Variable to Guessed copyValues(nx,ny); er = calcError(nx,ny); CPO(nx,ny); cout<<endl<<"Iteration # :"<<it<<"\tError:\t"<<er; if(it%AutoSave==0) { itoa(it,nm,10); ch=0; for(;nm[ch]!='\0';ch++){} nm[ch]='.'; nm[ch+1]='d'; nm[ch+2]='a'; nm[ch+3]='t'; nm[ch+4]='\0'; writetoFile(nx,ny,0,0,nm); cout<<"\nAutosave"<<endl; } it++; //system("cls"); }//END WHILE LOOP HERE itoa(it,nm,10); ch=0; for(;nm[ch]!='\0';ch++){} nm[ch]='.'; nm[ch+1]='d'; nm[ch+2]='a'; nm[ch+3]='t'; nm[ch+4]='\0'; end = clock(); // time reading for calculation on CPU time = (endstart_t)/double(CLK_TCK)*1000; cout<<"\nComputation Time = "<< time<<" ms\n"; writetoFile(nx,ny,time,er,nm); cout<<endl<<" DONE "<<endl; getchar(); } 

August 21, 2012, 06:00 

#32  
Senior Member
Join Date: Aug 2011
Posts: 271
Rep Power: 9 
Quote:
When Junaid started he was speaking about a twophase flow situation. As this kind of problem has a numerous potential sources of problems, I advised him to first check his NavierStokes solver in an incompressible classical problem like lidd driven cavity where it is easy to validate the code. Then he saw that he had some problems in his SIMPLE implementation. Now Junaid claimed that his code was runing well, but we don't know really where was the problem and how he overcomed it. But for me his linear system solver couldn't be the reason why it failed because, with the artificial compressibility algorithm it worked fine. Thus it means that the linear systeme solver worked fine to compute the velocity componnents. But as Arjun aptly mentionned it, his solver could run well with dirichlet boundary conditions, but not with neumann BC as it is required for pressure correction in the SIMPLE algorithm. It could be indeed the reason.. 

August 21, 2012, 07:30 

#33 
Member
Junaid Ahmad Khan
Join Date: Mar 2010
Location: Islamabad
Posts: 33
Rep Power: 9 
First let me clear that it is not fine but it gives similar results what LDC should give.
Problems: the residual of pressure is not going beyond 0.00012. it take urfp =0.0001 if i use 0.001 it diverges. at least 1st problem should not be there. So my next step is to remove that problem Thanks for all ur help i still need it though 

August 21, 2012, 08:22 

#34 
Senior Member
Join Date: Aug 2011
Posts: 271
Rep Power: 9 

August 21, 2012, 08:30 

#35  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 723
Rep Power: 19 
Quote:
if you have to go below 0.05 then either you have very very bad mesh or you are doing something really wrong. (by you i here mean a general person and not lefix) 

November 17, 2012, 14:36 
Dear Junaid  dividing by 2.0  bad programming style

#36 
Member
Michail
Join Date: Apr 2009
Location: Lithuania
Posts: 34
Rep Power: 10 
Dear Junaid  dividing by 2.0  bad programming style
Fw=A*(FP+FW)/2.0; Fe=A*(FE+FP)/2.0; Fs=A*(FSE+FSW)/2.0; Fn=A*(FNE+FNW)/2.0; cause usually multipication operation executed more rapidly than dividing 

May 12, 2013, 04:48 

#37 
New Member
hans
Join Date: Sep 2012
Posts: 7
Rep Power: 7 
Gents,
I'm running into the same problem. Looking at the reactions in this thread help me a bit, but still not there. If I look at eq. 6.36, 6.37 and the equations below that for the multiplier (d). I still end up with a problem. If I take the central differencing scheme from Chap. 5 for the central coefficient (pag 136): ap=aw+ae+FeFw; (Assume A and rho = 1 and D=0) aw=Fw/2; ae=Fe/2; Fw=uw;Fe=ue; ap=uw/2 ue/2 + uw/4 ue/4; When uw approaches ue (which happens if all works well in a 1D situation) ap goes to 0 and the multiplier (d) approaches infinity. Resulting in all sorts of problems. I'm not quite sure how the URF solves this issue. Can someone help me figuring out where I'm going wrong? 

May 12, 2013, 22:56 
same problem

#38 
New Member
David
Join Date: Feb 2013
Posts: 16
Rep Power: 6 
I am also implementing versteeg's book, the pressure correction got similar problem. The velocity distribution profile, pattern are all ok, but the magnitude value, if compared with slit flow theoretical value
u(y)=(dp/dx)H^2/8Miu*(1(y/0.5H)^2) the max velocity is always 50% higher than the theoretical number, I believe the problem is at pressure correction and don't know how to solve it. 

May 22, 2013, 04:14 

#40  
New Member
sina gilassi
Join Date: May 2013
Posts: 17
Rep Power: 6 
Quote:
based on the book "an introduction to fluid dynamic ..." page 142, there are two main equation of momentum, after I guess P(star), then the v and u (start) are calculated, my question is about the a.u(star) and the coefficient of these equation. how they are defined ?!!!! could you explain briefly thank you 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
SIMPLE algorithm in 3D cylindrical coordinates  zouchu  Main CFD Forum  1  January 20, 2014 18:02 
Finite Volume  SIMPLE Algorithm  Roger  Main CFD Forum  8  June 25, 2011 22:49 
SIMPLE algorithm confusion  lost.identity  Main CFD Forum  1  October 7, 2010 11:48 
SIMPLE OR SIMPLER algorithm  Sergio Costa  Main CFD Forum  2  July 29, 2007 06:44 
About Phase Coupled SIMPLE (PCSIMPLE) algorithm  Yan Kai  Main CFD Forum  0  April 18, 2007 03:48 