Implementation of Simple Algorithm
 

HI to All,

I am trying to solve Lid Driven Cavity problem using Finite Volume Method, for that i want to use Simple Algorithm. I have all ready solve LDC using Artificial Compressibility now i want to use SIMPLE Algorithm to solve velocity-pressure coupling, for that i am following Versteeg's book.

I am using steady navier stokes equations for LDC.

The procedure that i have adopted is given below:

1. Initialize u*,v*,p*
2. Solve Discretised momentum equation and Calculated u*,v*

NOW comes the problem. in step 3 when i try to calculate p'.
3.
the equation is
a[i,j]p'[i,j]=a[i-1,j]p'[i-1,j]+a[i+1,j]p'[i+1,j]+a[i,j-1]p'[i,j-1]+a[i,j+1]p'[i,j+1]+b'[i,j]

a[i-1,j]=(d A)[i-1,j]

a[i+1,j]=(d A)[i+1,j]

a[i,j-1]=(d A)[i,j-1]

a[i,j+1]=(d A)[i,j+1]

b'[i,j] = (u*A)[i-1,j]-(u*A)[i+1,j]+(v*A)[i,j+1]-(v*A)[i,j-1]
In this equation i need d to be calculated which is d=A/a where a is the centeral coefficient of velocity equation.
where i try to calculate this it some how a become zero which makes d=infinity. and p' become undefined or infinity. that cause problems in step 4 i.e.
4.
p[i,j]=p*[i,j]+p'[i,j]
u[i,j]=u*[i,j]+d[i,j]*(p'[i-1,j]-p'[i,j])
u[i,j]=u[i,j]+u[i,j]*(p'[i,j-1]-p'[i,j])
5. Set Boundary condition for u and v
6. Set Wall presure gradient to zero
7. Copy Values
p*=p
u*=u
v*=v
8. Check Convergence. and goto Step 1

My main problem is step 3 any sugesstions for that? (I am using c++)

Regards
Junaid

Hi Junaid

I do not know how you defined your notations but your implementation is a bit weird.

You should rather have something like:

ap[i,j]p'[i,j]=aw[i,j]p'[i-1,j]+ae[i,j]p'[i+1,j]+as[i,j]p'[i,j-1]+an[i,j]p'[i,j+1]+b'[i,j]

Yes your are correct yours make more sense. Mine your are same.
So any thoughts how to calculate d

Hi Junaid Ahmad!

Take a look here http://www.cfd-online.com/Wiki/Sampl...9_-_Fortran_90

Hope this will help

This is Pressure Correction File. First let me tell you that i have no working knowledge of FORTRAN So its hard for me to understand this. If I am not wrong then its your code. it will be very help full if you explain it and tell me which book you follow to develop your understanding. Or at least explain the Highlighted terms in this code.

Regards

Junaid

!************************************************* *********************
Subroutine Solve_Pressure_Correction
include 'icomm_1.f90'

DpU(2:NXmax,2:NYmax) = 1./Ap(2:NXmax,2:NYmax,1)
DpV(2:NXmax,2:NYmax) = 1./Ap(2:NXmax,2:NYmax,2)


!************************************************* ***************************************
! vertical faces East - e
Do 101 I=1,NXmax
Do 101 J=1,NYmax-1

If((i.ne.1).and.(i.ne.NXmax))then


DXc_e = Xc(i+1,j+1) - Xc(i,j+1)
S_e = Y(i ,j+1) - Y(i,j )

VOL_e = DXc_e * S_e

!-------------------------------------------------------------
APU_e = 0.5 * ( DpU(i,j+1) + DpU(i+1,j+1) )
Ul_e = 0.5 * ( F(i,j+1,1) + F(i+1,j+1,1) )
DPl_e = 0.5 * ( DPx_c(i,j+1) + DPx_c(i+1,j+1) )

DPx_e = ( F(i+1,j+1,4) - F(i,j+1,4) ) / DXc_e

U_e = Ul_e + APU_e * VOL_e * ( DPl_e - DPx_e)

!--------------------------------------------------------------
Con_e(i,j) = U_e * S_e


End If


Con_e(1,j) = 0.
Con_e(NXmax,j) = 0.

101 continue

!************************************************* ***************************************
! horisontal faces
Do 102 I=1,NXmax-1
Do 102 J=1,NYmax

If((j.ne.1).and.(j.ne.NYmax))then

DYc_n = Yc(i+1,j+1) - Yc(i+1,j)
S_n = X(i+1,j ) - X(i ,j)

VOL_n = DYc_n * S_n

!-----------------------------------------------------------
APV_n = 0.5 * ( DpV(i+1,j) + DpV(i+1,j+1) )
Vl_n = 0.5 * ( F(i+1,j,2) + F(i+1,j+1,2) )
DPl_n = 0.5 * ( DPy_c(i+1,j) + DPy_c(i+1,j+1) )


DPy_n = ( F(i+1,j+1,4) - F(i+1,j,4) ) / DYc_n

V_n = Vl_n + APV_n * VOL_n * ( DPl_n - DPy_n )
!-----------------------------------------------------------

Con_n(i,j) = V_n * S_n


End If


Con_n(i,1) = 0.
Con_n(i,NYmax) = 0.

102 continue

!************************************************* ***************************************
Summ = 0.

Do 200 I=2,NXmax
Do 200 J=2,NYmax

S_e = Y(i ,j) - Y(i,j-1)
S_n = X(i ,j) - X(i-1,j)


An(i,j) = 0.5 * ( DpV(i,j) + DpV(i,j+1) ) * S_n * S_n
As(i,j) = 0.5 * ( DpV(i,j) + DpV(i,j-1) ) * S_n * S_n
Ae(i,j) = 0.5 * ( DpU(i,j) + DpU(i+1,j) ) * S_e * S_e
Aw(i,j) = 0.5 * ( DpU(i,j) + DpU(i-1,j) ) * S_e * S_e

Ap(i,j,3) = (An(i,j) + As(i,j) + Ae(i,j) + Aw(i,j))

Sp(i,j,3) = -1. * ( Con_n(i-1,j) - Con_n(i-1,j-1) + Con_e(i,j-1) - Con_e(i-1,j-1) )

Summ = Summ + Sp(i,j,3) !
200 continue

write(*,*)'Sum',summ
!************************************************* **********************************
!************************************************* **********************************
write(*,*) Summ

write(*,*)'solve PP'
niter = 0


20 call Convergence_Criteria(3,Res_sum_before)
niter= niter + 1

Call TDMA_1(3)
call Convergence_Criteria(3,Res_sum_After)
If((abs(Res_sum_before-Res_sum_After).Ge.0.0000001).and.(niter.le.500).an d.(abs(Res_sum_After).ge.0.0001))go to 20
write(*,*)'Res_sum_before-Res_sum_After',Res_sum_before-Res_sum_After,niter




!************************************************* **********************************
!************************************************* **********************************
! velocities and pressure correction

Do 300 I=2,NXmax
Do 300 J=2,NYmax

DY = Y(i,j)-Y(i,j-1)
DX = X(i,j)-X(i-1,j)

PPe = 0.5 * ( F(i,j,3) + F(i+1,j,3) )
PPw = 0.5 * ( F(i,j,3) + F(i-1,j,3) )
PPn = 0.5 * ( F(i,j,3) + F(i,j+1,3) )
PPs = 0.5 * ( F(i,j,3) + F(i,j-1,3) )

F(i,j,1) = F(i,j,1) + (PPw - PPe) * DpU(i,j) * Dy
F(i,j,2) = F(i,j,2) + (PPs - PPn) * DpV(i,j) * Dx

F(i,j,4) = F(i,j,4) + 0.1 * F(i,j,3)

300 continue

!************************************************* **********************************
! correct convective fluxes through vertical East faces
Do 701 I=1,NXmax
Do 701 J=1,NYmax-1

If((i.ne.1).and.(i.ne.NXmax))then

Con_e(i,j) = Con_e(i,j) + Ae(i,j+1) * (F(i+1,j+1,3)-F(i,j+1,3) )

end if

if((i.eq.1).or.(i.eq.NXmax))Con_e(i,j)=0.

701 continue

! correct convective fluxes through horisontal North faces
Do 702 I=1,NXmax-1
Do 702 J=1,NYmax

If((j.ne.1).and.(j.ne.NYmax))then

Con_n(i,j) = Con_n(i,j) + An(i+1,j) * (F(i+1,j+1,3)-F(i+1,j,3) )

end if

if((j.eq.1).or.(j.eq.NYmax))Con_n(i,j)=0.

702 continue

!************************************************* **********************************
Do 501 I=1,NXmaxC

F(i,1,4) = F(i,2,4) !
F(i,NYmaxC,4) = F(i,NYmaxC-1,4)!


501 continue

Do 502 J=1,NYmaxC

F(1,j,4) = F(2,j,4) !
F(NXmaxC,j,4) = F(NXmaxC-1,j,4) !


502 continue

F(:,:,4) = F(:,:,4) - F(3,4,4)
!************************************************* **********************************



Return
End

Hi Junaid,

The code of Michail is really well writen and the notations are very understanable because intuitive. Greatings for him !!

You should first read a book on finite volume method and SIMPLE like the ones of Patankar, Peric or Versteek.
The one of Peric should be better because Michail used colocated variables arrangement based on Rhie and Chow interpolation. So for this purpose I think Peric's book is better oriented.

However just reading the code of Michail you understand what he is talking about...


Really get a book on FV , SIMPLE algorithm with Rhie and Chow interpolation and everything will be clear for you because the code is really clear!!

Thanks,

I am following verseek's book, and in it d(i-1,j)=A(i-1,j)/ap(i-1,j)

A(i-1,j) is area of west face
ap(i-1,j) is central coefficient of velocity
the problem is when i calculate ap by use the value of u* and v* and use it to find p' it some how blow up the velocity field, and when i use that to calculate ap which is ap = aw+ae+an+as+dF
For Hybrid scheme in that e.g. ae=Max(De-Fe/2,-Fe,0)
and Fe=0.5*(uP+uE)
De=Gama if dx=dy
so when velocity blowup everything went side ways.

So what should i use when i calculate ap do i use u* or u the corrected velocity.
Anyways i try to find this book of Peric and see how it works.

Dear Leflix

May I'll include Your explanations into the code? I guess it'l be better for future?

Regards
Michail

F(:,:,1) - U-velocity
F(:,:,2) - V-velocity
F(:,:,3) - pressure correction
F(:,:,4) - pressure

of course Michail, you are welcome!

SIMPLE FVM solution
 

I have a working Poiseulle flow problem based on Versteeg which matches nicely to the Hagen-Poiseulle velocity profile. It is extremey well annotated and culd be of great help.

It has under-relaxation and uses the JAcobi method, and was possible due to members of this forum such as Ztdep and houkensjtu.

The only problems that i have, are that the momentum resuduals do not get much better that 1e-3, and that the largest mesh that i am able to solve for (without divergence or stagnation) is about 34*34.... I would greatly appreciate any help, but for right now, you can look at my code...

Regards,
Michael

If you implemented versteeks method then would realy like to take a look at it. you may give me the like of send me the code at tojunaidahmad@hotmail.com

Regards
Junaid

Dear Junaid

I offer to take a look through M.Peric codes

ftp://ftp.springer.de/pub/technik/peric/

There are both staggered and collocated grids codes for SIMPLE-algorithm.

As concerned M.Peric book I offer to check torrents and previous posts on CFD-online.

Fight for CFD and CFD will fight for You :)

Best wishes
Michail

Thanks Michail, i need to know what is the best time to add under relaxation factor in the calculation of p, u, and v. now when i revisit the literature i find out that it is the under relaxation factor that cause divergence to me solution.

Regards
Junaid

were you using implicit under-relaxation for u and v? what values were you using for under-relaxation to make it diverge? I am having a similar issue...

what is the good method for under relaxation do i use






or use URF embedded in the discretised momentum equation

Regards
Junaid

If you look at Versteeg page 189, you will find that by applying the under-relaxation to the momentum equations yields the under-relaxed form of equations...6.36 and 6.37, they are not separate mehtods, and i believe that it is called implicit under-relaxation, and it the same form as that of peric.

AS for pressure under-relaxation, I understand it such that only the correction is under-relaxed... i.e. pnew = p* +alphap*p'

I don't understand i did exactly what is written in book. it didn't work. I did the same think using artificial compressibility and it worked but using SIMPLE method it did not. the only thing left for me is to paste the code a let you guys decide what is wrong with it.

But first i have to go through one more time.

It sounds that easy! ;)

Did you get the code that i sent? I have now upgraded it to SIMPLEC also, and it is working well.

Yes i got the code . and its approach and mine is the same i.e. hi-bride scheme.


why we need SIMPLE Algorithm when We have Artificial compressibility?:confused:

thank you would be nice.


as for artificial compressibility, that is outside of my knowledge.

Artificial compressibility is designed and works only for steady problems.

if your artificial compressiblity method works and SIMPLE does not work then i would suggest have a look at your linear system solver, that probably is not working well for SIMPLE case. Thats because the nature of linear system in both cases is very different.

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.

i think this part is never ruled out.

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

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.

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.

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.

As I mentioned it previously there is no linear system to solve for pressure in the artificial compressibility method.

That's why I think it is more stable, pressure was calculate explicitly.
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.

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

//Under-Relaxation 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};
//---------------------------------

//U-velocity, V-velocity, 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 - sub-diagonal (means it is the diagonal below the main diagonal) -- indexed from 1..n-1
* b - the main diagonal
* c - sup-diagonal (means it is the diagonal above the main diagonal) -- indexed from 0..n-2
* v - right part
* x - the answer
*/
for (int i = 1; i < n; i++)
{
double m = a[i]/b[i-1];
b[i] = b[i] - m*c[i-1];
v[i] = v[i] - m*v[i-1];
}

x[n-1] = v[n-1]/b[n-1];

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[Nx-1][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][Nx-1];//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[col-1][row]*A;
aE=du[col][row]*A;
aS=dv[col][row]*A;
aN=dv[col][row-1]*A;


aP=aW+aE+aS+aN;

bpw=us[col-1][row]*A;
bpe=us[col][row]*A;
bps=vs[col][row]*A;
bpn=vs[col][row-1]*A;

bp=bpw-bpe+bps-bpn;


pp[col][row]=(
aW*pp[col-1][row]
+ aE*pp[col+1][row]
+ aS*pp[col][row+1]
+ aN*pp[col][row-1]
+ bp
)/aP;

//a(I,J) p'(I,J)=a(I+1,J) x p'(I+1,J ) + a(I-1,J) x p'(I-1,J ) + a(I,J+1) x p'(I,J+1) + a(I,J-1) x p'(I,J-1) + 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)
{
//U-Velocity Corrector
//----------------------------------

float d=0.0;
for(int row=1; row<Ny; row++)
{
for(int col=1; col<Nx-1; 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))
}
}
}

//U-under relaxation
void U_urf(int Nx,int Ny,float A,float G)
{
//U-Velocity 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==0||col==(Nx-1))
{
a1[sb]=0.0;
c1[sb]=0.0;
b1[sb]=1.0;
v1[sb]=u[col][row];
}
else if(col>0&&col<(Nx-1))
{
FP=u[col][row];
FW=u[col-1][row];
FE=u[col+1][row];
FSE=v[col+1][row];
FSW=v[col][row];
FNE=v[col+1][row-1];
FNW=v[col][row-1];

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=Fe-Fw+Fn-Fs;

aw=MAX(G+Fw/2.0,MAX(Fw,0.0));
ae=MAX(G-Fe/2.0,MAX(-Fe,0.0));
as=MAX(G+Fs/2.0,MAX(Fs,0.0));
an=MAX(G-Fn/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][row-1]
+ A*(p[col][row]-p[col+1][row])
+ ((1-urfu)*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 +[(1-alpha_u)a(i,j)/alpha_u] x u^(n-1)_(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)
{
//V-Velocity Corrector
float v_calc=0.0,d=0.0;
for(int col=1; col<Nx; col++)
{
for(int row=1; row<Ny-1; 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))
}
}
}


//V-under 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==0||row==(Ny-1))
{
a2[sb]=0.0;
c2[sb]=0.0;
b2[sb]=1.0;
v2[sb]=v[col][row];
}
else if(row>0&&row<(Ny-1))
{
FSW=(u[col-1][row+1]);
FNW=(u[col-1][row]);
FNE=(u[col][row]);
FSE=(u[col][row+1]);
FP=(v[col][row]);
FS=(v[col][row+1]);
FN=(v[col][row-1]);

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=Fe-Fw+Fn-Fs;

D=G;


aw=MAX(D+Fw/2.0,MAX(Fw,0.0));
ae=MAX(D-Fe/2.0,MAX(-Fe,0.0));
as=MAX(D+Fs/2.0,MAX(Fs,0.0));
an=MAX(D-Fn/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[col-1][row]
+ ae*v[col+1][row]
+ A*(p[col][row+1]-p[col][row])
+ ((1-urfv)*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 +[(1-alpha_v)a(I,j)/alpha_v] x v^(n-1)_(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][Ny-1];

}

for (int row=0; row<=Ny; row++)
{
//West
p[0][row] = p[1][row];

//East
p[Nx][row] = p[Nx-1][row];
}

}


void calcUs(int Nx,int Ny,float A,float G)
{
/***********Calculate u-Velocity 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==0||col==(Nx-1))
{
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<(Nx-1))
{
FP=us[col][row];
FW=us[col-1][row];
FE=us[col+1][row];
FSE=vs[col+1][row];
FSW=vs[col][row];
FNE=vs[col+1][row-1];
FNW=vs[col][row-1];

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=Fe-Fw+Fn-Fs;

aw=MAX(G+Fw/2.0,MAX(Fw,0.0));
ae=MAX(G-Fe/2.0,MAX(-Fe,0.0));
as=MAX(G+Fs/2.0,MAX(Fs,0.0));
an=MAX(G-Fn/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][row-1]
+ 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 v-Velocity 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==0||row==(Ny-1))
{
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<(Ny-1))
{
FSW=(us[col-1][row+1]);
FNW=(us[col-1][row]);
FNE=(us[col][row]);
FSE=(us[col][row+1]);
FP=(vs[col][row]);
FS=(vs[col][row+1]);
FN=(vs[col][row-1]);

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=Fe-Fw+Fn-Fs;

D=G;


aw=MAX(D+Fw/2.0,MAX(Fw,0.0));
ae=MAX(D-Fe/2.0,MAX(-Fe,0.0));
as=MAX(D+Fs/2.0,MAX(Fs,0.0));
an=MAX(D-Fn/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[col-1][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/(nx-1);

float deltaY = Ly/(ny-1);

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 = (end-start_t)/double(CLK_TCK)*1000;
cout<<"\nComputation Time = "<< time<<" ms\n";

writetoFile(nx,ny,time,er,nm);
cout<<endl<<" DONE "<<endl;

getchar();
}

When Junaid started he was speaking about a two-phase flow situation.
As this kind of problem has a numerous potential sources of problems, I advised him to first check his Navier-Stokes 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..

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

Using a URF_p <0.1-0.2 is very pathologic
It's not acceptable especially for lidd driven cavity problem.

+1


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)

Dear Junaid - dividing by 2.0 - bad programming style
 

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

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+Fe-Fw; (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? :confused:

same problem
 

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.

northfly,

Do you run into the same problem with the multiplier (d) that approaches infinity when all neighbour cells (velocity )have the same value? How have you solved this?

Kind regards

Dear ahmad,

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