CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

continuity equation satisfaction

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

Reply
 
LinkBack Thread Tools Display Modes
Old   July 17, 2010, 14:06
Default continuity equation satisfaction
  #1
New Member
 
sayyed sadegh forghani dehnavi
Join Date: Jul 2010
Location: tehran
Posts: 2
Rep Power: 0
sadeghforghani is on a distinguished road
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!!! 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;
}

Last edited by sadeghforghani; July 17, 2010 at 14:25.
sadeghforghani is offline   Reply With Quote

Old   February 4, 2011, 10:56
Default
  #2
Senior Member
 
maysmech's Avatar
 
Join Date: Jan 2010
Posts: 347
Blog Entries: 2
Rep Power: 8
maysmech is on a distinguished road
Quote:
Originally Posted by sadeghforghani View Post
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!!! 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;
}
Hi Sadegh,
Is your problem solved?
maysmech is offline   Reply With Quote

Old   February 5, 2011, 04:00
Default
  #3
New Member
 
sayyed sadegh forghani dehnavi
Join Date: Jul 2010
Location: tehran
Posts: 2
Rep Power: 0
sadeghforghani is on a distinguished road
Quote:
Originally Posted by maysmech View Post
Hi Sadegh,
Is your problem solved?
Dear maymech,

Yes, I found error!

Thanks,
sadeghforghani is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Calculation of the Governing Equations Mihail CFX 7 September 7, 2014 06:27
How to write k and epsilon before the abnormal end xiuying OpenFOAM Running, Solving & CFD 8 August 27, 2013 15:33
IcoFoam parallel woes msrinath80 OpenFOAM Running, Solving & CFD 9 July 22, 2007 02:58
continuity equation Rafal Main CFD Forum 4 November 29, 2006 10:27
Could anybody help me see this error and give help liugx212 OpenFOAM Running, Solving & CFD 3 January 4, 2006 19:07


All times are GMT -4. The time now is 00:04.