CFD Online Discussion Forums - continuity equation satisfaction

Hi:)

I have a c++ code(you can see it below) to solve N.S axisymmetric equations in a pipe with staggered lattice. but after the convergency I found that there was not a continuity satisfaction!!! :mad:because the inlet flow flux was more than the outlet.
please help me correcting it.

regards,

code:
#include<iostream.h>
#include<stdio.h>
#include<math.h>
#include<fstream.h>

long double un[500],u[500][500],ustar[500][500],unew[500][500],dur[500][500],duz[500][500],d2ur[500][500],d2uz[500][500];
long double vn[500],v[500][500],vstar[500][500],vnew[500][500],dvr[500][500],dvz[500][500],d2vr[500][500],d2vz[500][500];
long double wn[500],w[500][500],wstar[500][500],wnew[500][500],dwr[500][500],dwz[500][500],d2wr[500][500],d2wz[500][500];

long double uconv[500][500],udiff[500][500],erru[500][500], p[500][500], miu[500][500], ro[500][500], errp[500][500];
long double vconv[500][500],vdiff[500][500],errv[500][500],dpr[500][500],d2pr[500][500], nu ,dustarr[500][500];
long double wconv[500][500],wdiff[500][500],errw[500][500],dpz[500][500],d2pz[500][500],ronew[500][500],dwstarz[500][500];

long double gridc[500][2],eru,erp,pnew,gr,r ;
long double grid[500][500][2],erv,roa,miua,gz,rx,x,xc,delr;
long double delta[500][500][2],erw,row,miuw,gt,ry,y,yc,delz;

long double a1,a4,a5,a6,a7,a8,a9,a10;
long double p1,p2,p3,p4,coeffp;
long double delteta,d,dmax,delt,teta,error,maxerror,gradient,p i,test2;

int i,j,k,l,m,n,levelset;
int gridnod[500][2],a2,a3,iter,piter,t,test,nc,nodc,maxiter;
int slnh1,slnh2,slnv1,slnv2;
/************************************************** *******************/
void constants();
void gridd();
void griddc();
void initial();
void settingboundaries();
/************************************************** ********************/
int main()
{
constants();
gridd();
griddc();
initial();
iter=0;
error=1.0;
cout<<"dx="<<delr<<" dy="<<delz<<endl;
while(error>maxerror)
{
settingboundaries();

ofstream out;
out.open("out.dat");
out<<"title= output"<<endl;
out<<"variables= x , y , u , w ,p"<<endl;
out<<"zone i="<<n<<" j="<<m<<endl;
for(i=3;i<=(2*m+1);i=i+2)
for(j=3;j<=(2*n+1);j=j+2)
out<<grid[i][j][1]<<" "<<grid[i][j][2]<<" "<<(u[i][j-1]+u[i][j+1])/2.0<<" "<<(w[i-1][j]+w[i+1][j])/2.0<<" "<<(p[i-1][j-1]+p[i-1][j+1]+p[i+1][j-1]+p[i+1][j+1])/4.0<<endl;
out.close();
iter++;
for(i=5;i<=(2*m-1);i=i+2)//uuuuuu
for(j=4;j<=(2*n);j=j+2)
{
dur[i][j]=(u[i+2][j]-u[i-2][j])/(4*delr);
duz[i][j]=(u[i][j+2]-u[i][j-2])/(4*delz);
d2ur[i][j]=(u[i+2][j]-2*u[i][j]+u[i-2][j])/(4*delr*delr);
d2uz[i][j]=(u[i][j+2]-2*u[i][j]+u[i][j-2])/(4*delz*delz);
}
for(i=5;i<=(2*m-1);i=i+2)//uuuuuu
for(j=4;j<=(2*n);j=j+2)
{
r=((2*m+1)-i)*delr;//check shavad
uconv[i][j]=-u[i][j]*dur[i][j]-/*w[i][j]*/((w[i-1][j-1]+w[i-1][j+1]+w[i+1][j-1]+w[i+1][j+1])/4.0)*duz[i][j];
udiff[i][j]=d2ur[i][j]+dur[i][j]/r-u[i][j]/(r*r)+d2uz[i][j];
}

for(i=4;i<=(2*m);i=i+2)//wwwwwww
for(j=5;j<=(2*n-1);j=j+2)
{
dwr[i][j]=(w[i+2][j]-w[i-2][j])/(4*delr);
dwz[i][j]=(w[i][j+2]-w[i][j-2])/(4*delz);
d2wr[i][j]=(w[i+2][j]-2*w[i][j]+w[i-2][j])/(4*delr*delr);
d2wz[i][j]=(w[i][j+2]-2*w[i][j]+w[i][j-2])/(4*delz*delz);
}
for(i=4;i<=(2*m);i=i+2)//wwwwwww
for(j=5;j<=(2*n-1);j=j+2)
{
r=((2*m+1)-i)*delr;//check shavad
wconv[i][j]=-/*u[i][j]*/((u[i-1][j-1]+u[i-1][j+1]+u[i+1][j-1]+u[i+1][j+1])/4.0)*dwr[i][j]-w[i][j]*dwz[i][j];
wdiff[i][j]=d2wr[i][j]+dwr[i][j]/r+d2wz[i][j];
}

for(i=5;i<=(2*m-1);i=i+2)//uuuuuu
for(j=4;j<=(2*n);j=j+2)
ustar[i][j]=delt*(uconv[i][j]+nu*udiff[i][j]+gr)+u[i][j];

for(i=4;i<=(2*m);i=i+2)//wwwwwww
for(j=5;j<=(2*n-1);j=j+2)
wstar[i][j]=delt*(wconv[i][j]+nu*wdiff[i][j]+gz)+w[i][j];
// printf("uconv=%5.5f , vconv=%5.5f , wconv=%5.5f , udiff=%5.5f , vdiff=%5.5f , wdiff=%5.5f , ustar[%d][%d]=%5.5f , vstar[%d][%d]=%5.5f , wstar[%d][%d]=%5.5f\n",uconv,vconv,wconv,udiff,vdiff,wdiff,i,j, ustar[i][j],i,j,vstar[i][j],i,j,wstar[i][j]);
/************************************************** ***********************/
piter=1;
erp=1.0;
while((erp>0.000000000000001)&&(piter<2000000))
{
erp=0.0;
for(i=2;i<=(2*m+2);i=i+2)
{
p[i][2]=p[i][4];
p[i][2*n+2]=p[i][2*n];
}
for(i=2;i<=(2*m+2);i=i+2)
{
wstar[i][1]=w[i][1];
wstar[i][3]=w[i][3];
wstar[i][2*n+1]=w[i][2*n+1];
wstar[i][2*n+3]=w[i][2*n+3];
}
for(i=1;i<=(2*m+3);i=i+2)
{
ustar[i][2]=u[i][2];
ustar[i][2*n+2]=u[i][2*n+2];
}

for(j=2;j<=(2*n+2);j=j+2)
{
p[2][j]=p[4][j];
p[2*m+2][j]=p[2*m][j];
}
for(j=2;j<=(2*n+2);j=j+2)
{
ustar[1][j]=u[1][j];
ustar[3][j]=u[3][j];
ustar[2*m+1][j]=u[2*m+1][j];
ustar[2*m+3][j]=u[2*m+3][j];
}
for(j=1;j<=(2*n+3);j=j+2)
{
wstar[2][j]=w[2][j];
wstar[2*m+2][j]=w[m][j];
}
for(i=2;i<=(2*m+2);i=i+2)
for(j=2;j<=(2*n+2);j=j+2)
errp[i][j]=0;

/* for(i=5;i<=(2*m-1);i=i+2)//uuuuuu
for(j=4;j<=(2*n);j=j+2)
{
r=((2*m+1)-i)*delr;//check shavad
dustarr[i][j]=(ustar[i+2][j]-ustar[i-2][j])/(4*delr);
}
for(i=4;i<=(2*m);i=i+2)//wwwwwww
for(j=5;j<=(2*n-1);j=j+2)
{
r=((2*m+1)-i)*delr;//check shavad
dwstarz[i][j]=(wstar[i][j+2]-wstar[i][j-2])/(4*delz);
}
*/ for(i=4;i<=(2*m);i=i+2)//pppppppppppppppppp
for(j=4;j<(2*n);j=j+2)
{
r=((2*m+1)-i)*delr;//check shavad
coeffp=-2.0/(4*delr*delr)-2.0/(4*delz*delz);//hameye delha zarb dar 2 shodand
p1=/*ustar[i][j]*/((ustar[i-1][j]+ustar[i+1][j])/2.0)/r+
/*dustarr[i][j]*/((ustar[i+1][j]-ustar[i-1][j])/(2.0*delr))+
/*dwstarz[i][j]*/((wstar[i][j+1]-wstar[i][j-1])/(2.0*delz));
p2=(p[i+2][j]-p[i-2][j])/(4*delr*r);
p3=(p[i+2][j]+p[i-2][j])/(4*delr*delr);
p4=(p[i][j+2]+p[i][j-2])/(4*delz*delz);
pnew=((row/delt)*p1-p2-p3-p4)/(coeffp);
errp[i][j]=fabs((pnew-p[i][j])/pnew);
erp=erp+errp[i][j]*errp[i][j];
p[i][j]=pnew;
}
erp=sqrt(erp);
// if(piter%100000==0)
// printf("piter=%10.10d , erp=%20.17f\n",piter,erp);
piter++;
}
// if(iter%10==0)
// printf("piter=%10.10d , erp=%20.17f\n",piter,erp);
// if(iter%100==0)
// printf("p @ %d converged & piter=%d & erp=%20.17f\n",iter,piter,erp);
/************************************************** *************************************/
eru=0.0;
erv=0.0;
erw=0.0;
for(i=5;i<=(2*m-1);i=i+2)//uuuuuu
for(j=4;j<=(2*n);j=j+2)
erru[i][j]=0.0;
for(i=4;i<=(2*m);i=i+2)//wwwwwww
for(j=5;j<=(2*n-1);j=j+2)
errw[i][j]=0.0;

/* for(i=4;i<=(2*m);i=i+2)//ppppppp
for(j=4;j<(2*n);j=j+2)
{
dpr[i][j]=(p[i+2][j]-p[i-2][j])/(4*delr);
dpz[i][j]=(p[i][j+2]-p[i][j-2])/(4*delz);
}
*/ test2=0.0;
for(i=5;i<=(2*m-1);i=i+2)//uuuuuu
for(j=4;j<=(2*n);j=j+2)
{
/* if(ustar[i][j]>test2)
{
test2=ustar[i][j];
cout<<"test2="<<test2<<endl;
}
*/ unew[i][j]=ustar[i][j]-(delt/row)*((p[i+1][j]-p[i-1][j])/(2.0*delr));/*dpr[i][j]*/
erru[i][j]=fabs(u[i][j]-unew[i][j])/delt;
u[i][j]=unew[i][j];
eru=eru+erru[i][j]*erru[i][j];
}

for(i=4;i<=(2*m);i=i+2)//wwwwwww
for(j=5;j<=(2*n-1);j=j+2)
{
wnew[i][j]=wstar[i][j]-(delt/row)*((p[i][j+1]-p[i][j-1])/(2.0*delz));/*dpz[i][j]*/
errw[i][j]=fabs(w[i][j]-wnew[i][j])/delt;
w[i][j]=wnew[i][j];
erw=erw+errw[i][j]*errw[i][j];

}

// cout<<"eru="<<eru<<" erw="<<erw<<endl;
error=sqrt(eru+erv+erw);
if(iter%100==0)
printf("iterate number=%5.5d error=%20.17f\n",iter,error);
//error=1;
}
cout<<"final iter number="<<iter<<" , last error="<<error<<endl;
cout<<"the end!"<<endl;
return 0;
}

/************************************************** ***********************************/
//edited
void settingboundaries()
{
for(i=1;i<=(2*m+3);i=i+2)
{
u[i][2]=-u[i][4];
u[i][2*n+2]=u[i][2*n];
}
for(i=2;i<=(2*m+2);i=i+2)
{
w[i][3]=0.01;
w[i][1]=0.01;
w[i][2*n+1]=w[i][2*n-1];
w[i][2*n+3]=w[i][2*n-1];
}
for(j=2;j<=(2*n+2);j=j+2)
{
u[3][j]=0.0;
u[1][j]=0.0;
u[2*m+1][j]=0.0;
u[2*m+3][j]=0.0;
}
for(j=1;j<=(2*n+3);j=j+2)
{
w[2][j]=-w[4][j];
w[2*m+2][j]=w[2*m][j];
}
return;
}
/************************************************** ***********************/
//edited
void gridd()
{
for(j=3;j<=(2*n+1);j++)
for(i=3;i<=(2*m+1);i++)
{
grid[i][j][1]=(i-3)*delr-(x/2.0);
grid[i][j][2]=(j-3)*delz-(y/2.0);
}
return;
}
//************************************************** **********
//edit nadasht
void griddc()
{
delteta=pi/(nc-1);
teta=pi/2.0;
for(i=1;i<=nc;i++)
{
gridc[i][1]=(d/2.0)*cos(teta);
if(gridc[i][1]>0)
gridc[i][1]=0.0;
gridc[i][2]=(d/2.0)*sin(teta)-(y/4.0);
teta=teta+delteta;
}
return;
}
//################################################## ########
//edited
void initial()
{
for(i=1;i<=(2*m+3);i=i+2)
for(j=2;j<=(2*n+2);j=j+2)
{
u[i][j]=0.0;
ustar[i][j]=0.0;
}
for(i=2;i<=(2*m+2);i=i+2)
for(j=1;j<=(2*n+3);j=j+2)
{
w[i][j]=0.0;
wstar[i][j]=0.0;
}
for(i=2;i<=(2*m+2);i=i+2)
for(j=2;j<=(2*n+2);j=j+2)
p[i][j]=10000;
return;
}
//************************************************** ********************/
void constants()
{
maxiter=1000;
maxerror=0.000000000001;
m=20;//tedade nodha
n=2*m;
dmax=0.008;
d=0.5*dmax;
levelset=0;//0 yani initial va 1 yani az edame hale ghabli
pi=4.0*atan(1.0);
nc=60;
x=5.0*dmax;
y=10.0*dmax;
delr=x/(2*(2*(m-1)));
delz=y/(2*(n-1));
gr=0.0;
gt=0.0;
gz=0.0;
delt=0.01;
miuw=0.001002;
row=1000.0;
nu=miuw/row;
return;
}

Quote:

Originally Posted by sadeghforghani (Post 267827)

Hi Sadegh,
Is your problem solved?