sci.physics.computational.fluid-dynamics

As in my last post i am trying to develop a small code to simulate
turbulent flow in pipe using k-eps model
I have added a new post because i am now changed my approach i am
solving the vel components and k and eps equations explicitly and
pressure term implicitly
i tried diifferent things for boundary conditions but still not
working one more problem which i faced was i am checking for
continuity eq and iterate till the max value reaches almost zero but
the problem is that all of them never converge in fact the values
increase after some time
what should i do ...adding my code for your ref......
#include
#include
//#include
#include
#include
#define i0 100
#define j0 10
#define pos0 5
#define t0 10
# define PI 3.1429
#define out 0

double cont();
//double tau(double,double,double,double);
double funE(int,int,int,int);
double funF(int,int,int,int);
double funG(int,int,int,int);
double implicitpart(double);
void keps();
void explicitpart(int);
int i,j,pos,n,iter;
double
rho,ptime,visc,eij,cu,r,dt,dx,dr,sigmak,sigmaeps,c1eps,c2eps,ut,rad,r1,dt1,beta,dmax;
double du[pos0],du1[pos0],u[i0][j0][t0],v[i0][j0][t0],w[i0][j0][t0];
double k[i0][j0][t0],eps[i0][j0][t0],pr[i0][j0][t0],visct[i0][j0]
[t0],d[i0][j0],ubar[i0][j0][t0],vbar[i0][j0][t0];
int main()
{
double dmax1,dtest,x=0,flag;
ofstream file("file.xls");
ofstream dma("dmax.xls");
//constants for the flow.....considering dx=dr for grids
rad=0.2;
r=0.5;
dx=dr=0.02;
dt=dr*r;
beta=1/(2*dt*(1/(dx*dx)+1/(dr*dr)));
dt1=dt;
rho=1;
visc=2e-5;
//constants for k-eps model
cu=0.09;
c1eps=1.44;
c2eps=1.92;
sigmak=1.0;
sigmaeps=1.0;

//Initializing the varialbles for t=0
n=0;
for(i=0;i {
for(j=0;j {
u[i][j][n]=1.0;
v[i][j][n]=0.0;
pr[i][j][n]=101325.0;
k[i][j][n]=(1/10.0)*(u[i][j][0]*u[i][j][0]);
eps[i][j][n]=pow(k[i][j][0],1.5)/dr;
visct[i][j][n] = cu*(pow(k[i][j][n],1.5))/eps[i][j][n]
+visc;
}
}
dtest=0.005*u[0][0][0]/dx;

//loop for calculation
for(n=0;n {
flag=0;
file< dmax1=1.0;
iter=0;
while((fabs(dmax1))>= dtest)
{
keps();//k-eps equation
if(flag==0)
{
explicitpart(n);//momemtum equation without pressure term and
du/dt
flag=1;
}
dmax1=implicitpart(beta);//eveluating pressure term implicitly
and then calculating u using that pressure
iter++;
if(iter>100)//just to avoid excess computation and break loop in
that case
break;
}
for(i=0;i {
for(j=0;j {
file<[n]<<"\t"< }
}
}
//display

file.close();
dma.close();
return 0;
}

double funE(int i1, int j1, int pos1, int t1)
{
if(pos1==0)
return (u[i1][j1][t1]*u[i1][j1][t1]-(visct[i1][j1][t1]/dr)*(u[i1+1][j1]
[t1]-u[i1][j1][t1]));//evaluating dE/dx for x dir mom eq

if(pos1==1)
return (u[i1][j1][t1]*v[i1][j1][t1]-(visct[i1][j1][t1]/dr)*(v[i1+1][j1]
[t1]-v[i1][j1][t1])); //evaluating dE/dx for r dir mom eq

if(pos1==2)
return (rho*k[i1][j1][t1]*u[i1][j1][t1]-(visct[i1][j1][t1]/
sigmak)*(k[i1+1][j1][t1]-k[i1][j1][t1])/dr);//evaluating dE/dx for k
eq

if(pos1==3)
return (rho*eps[i1][j1][t1]*u[i1][j1][t1]-(visct[i1][j1][t1]/
sigmaeps)*(eps[i1+1][j1][t1]-eps[i1][j1][t1])/dr);//evaluating dE/dx
for eps eq

}

double funF(int i2,int j2, int pos2, int t2)
{
if(pos2==0)
return (u[i2][j2][t2]*v[i2][j2][t2]-(visct[i2][j2][t2]/dr*(u[i2][j2+1]
[t2]-u[i2][j2][t2]))-(visct[i2][j2][t2]*u[i2][j2][t2]/r1));//
evaluating dF/dr for x dir mom eq

if(pos2==1)
return (v[i2][j2][t2]*v[i2][j2][t2]-(visct[i2][j2][t2]/dr*(v[i2][j2+1]
[t2]-v[i2][j2][t2]))-(visct[i2][j2][t2]*v[i2][j2][t2]/r1));//
evaluating dF/dr for r dir mom eq

if(pos2==2)
return (rho*k[i2][j2][t2]*u[i2][j2][t2]-(visct[i2][j2][t2]/
sigmak)*(k[i2][j2+1][t2]-k[i2][j2][t2])/dr);//evaluating dF/dr for k
eq

if(pos2==3)
return (rho*eps[i2][j2][t2]*u[i2][j2][t2]-(visct[i2][j2][t2]/
sigmaeps)*(eps[i2][j2+1][t2]-eps[i2][j2][t2])/dr);//evaluating dF/dr
for eps eq

}

double funG(int i3,int j3, int pos3, int t3)
{
double mixed;
mixed=(1/(dr*dr))*((u[i3][j3+1][t3]-u[i3][j3][t3]+v[i3+1][j3][t3]-v[i3]
[j3][t3])*(u[i3][j3+1][t3]-u[i3][j3][t3]+v[i3+1][j3][t3]-v[i3][j3]
[t3])/4);
eij=(1/(dr*dr))*((2*(u[i3+1][j3][t3]-u[i3][j3][t3])*(u[i3+1][j3][t3]-
u[i3][j3][t3]))+(2*(v[i3][j3+1][t3]-v[i3][j3][t3])*(v[i3][j3+1][t3]-
v[i3][j3][t3])))+ mixed;
//evaluating the productiong term in k-eps eq
if(pos3==0)
return (visc*u[i3][j3][t3]/(r1*r1));//evaluating G for x dir mom eq

if(pos3==1)
return (visc*v[i3][j3][t3]/(r1*r1));//evaluating G for r dir mom eq

if(pos3==2)
return (2*visct[i3][j3][t3]*eij*eij-rho*eps[i3][j3][t3]);//evaluating
G for k eq

if(pos3==3)
return (c1eps*eps[i3][j3][t3]*2*visct[i3][j3][t3]*eij*eij-
c2eps*rho*eps[i3][j3][t3]*eps[i3][j3][t3]/k[i3][j3][t3]);//evaluating
G for eps eq
}

double implicitpart(double beta)
{int a,b;
for(i=1;i {
r1=dr*2;
for(j=1;j {
if(iter>0)
{
u[i][j][n]=u[i][j][n+1];
v[i][j][n]=v[i][j][n+1];
u[i-1][j][n]=u[i-1][j][n+1];
v[i][j-1][n]=v[i][j-1][n+1];
}
d[i][j]=((u[i][j][n]-u[i-1][j][n])/dr+(v[i][j][n]-v[i][j-1][n])/
dr+(v[i][j][n]/r1));//continuity equation
if ((i<1)&& ((fabs(d[i][j]))>dmax))
dmax = d[i][j];//max value of continuity eq to compare for
iteration
//if(n>0&&j==8) {cout<"<[j]< r1=r1+dr;
}
}
for(i=1;i {
for(j=1;j {
pr[i][j][n]=pr[i][j][n]-beta*d[i][j];//pressure computation
// if((i==1)&&(j==8)&&(n>0)){cout<"< }
}

for(i=0;i {
r1=dr;
for(j=0;j {
if((i>0)&&(j>0))
{
u[i][j][n+1] = ubar[i][j][n]+(dt/dx)*(pr[i][j][n]-pr[i-1][j]
[n]);//velocity computations

v[i][j][n+1] = vbar[i][j][n]+(dt/dr)*(pr[i][j][n]-pr[i][j-1]
[n]);
}
//boundaries for vel compoents
if((i==2)&&(j!=j0-1))
{//inlet boundary conditions
u[0][j][n+1] = fabs(sin(PI*0.1*n));
v[0][j][n+1] = 2*u[1][j][n+1]-u[2][j][n+1];
pr[0][j][n+1] = 2*pr[1][j][n+1]-pr[2][j][n+1];
}
if(j==1)
{//axial boundary case
u[i][j-1][n+1]=u[i][j][n+1];
v[i][j-1][n+1]=v[i][j][n+1];
}
if((i>i0-3)&&(j!=j0-1))
{//outlet boundary conditions
u[i][j][n+1]=2*u[i-2][j][n+1]-u[i-3][j][n+1];
v[i][j][n+1]=2*v[i-2][j][n+1]-v[i-3][j][n+1];
}
if((j==j0-1)&&(i>0)&&(i {//wall boundary conditions
u[i][j][n+1]=0.0;
v[i][j][n+1]=0.0;
visct[i][j][n+1] = 0.0;
ut=sqrt(visct[i][j][n+1]/rho*(u[i][j0-2][n+1]-u[i][j0-1][n
+1])/dr);
k[i][j0-1][n+1]=0.0;
eps[i][j0-1][n+1]=0.0;
}
r1=r1+dr;
}
}

return dmax;//for comparsion
}

void explicitpart(int n)
{
for(i=1;i {
r1=dr;
for(j=0;j {
du[0] = u[i][j][n];
du[1] = v[i][j][n];
ubar[i][j][n] = du[0] - (r*(funE(i+1,j,0,n) - funE(i,j,0,n))) -
(r*(funF(i,j+1,0,n) - funF(i,j,0,n))) + (dt*funG(i,j,0,n));//explicit
part of
// if(j==8&&n>0){cout<"< vbar[i][j][n] = du[1] - (r*(funE(i+1,j,1,n) - funE(i,j,1,n))) -
(r*(funF(i,j+1,1,n) - funF(i,j,1,n))) + (dt*funG(i,j,1,n));//vel
compoennts at n+1
r1=r1+dr;
}
}
}

void keps()
{
for(i=0;i {
r1=dr;
for(j=0;j {
du[2] = k[i][j][n]*rho;
du[3] = eps[i][j][n]*rho;

du1[2] = du[2] - (r*(funE(i+1,j,2,n) - funE(i,j,2,n))) -
(r*(funF(i,j+1,2,n) - funF(i,j,2,n))) + (dt*funG(i,j,2,n));//
evaluation of k at n+1
du1[3] = du[3] - (r*(funE(i+1,j,3,n) - funE(i,j,3,n))) -
(r*(funF(i,j+1,3,n) - funF(i,j,3,n))) + (dt*funG(i,j,3,n));//
evaluation of eps at n+1

k[i][j][n+1] = du1[2]/rho;
eps[i][j][n+1] = du1[3]/rho;

r1=r1+dr;

//boundaries for k eps eq
if((i==2)&&(j!=j0-1))
{//inlet boundary conditions
k[0][j][n+1] = 2*u[1][j][n+1]-u[1][j][n+1];
eps[i-2][j][n+1]=pow(k[0][j][n+1],1.5);
}
if(j==1)
{//inlet axial conditions
k[i][j-1][n+1]=k[i][j][n+1];
eps[i][j-1][n+1]=eps[i][j][n+1];
}
if((i>i0-3)&&(j!=j0-1))
{//outlet boundary conditions
k[i][j][n+1]=2*k[i-2][j][n+1]-k[i-3][j][n+1];
eps[i][j][n+1]=2*eps[i-2][j][n+1]-eps[i-3][j][n+1];
}
visct[i][j][n+1] = cu*(pow(k[i][j][n+1],1.5))/eps[i][j][n
+1]+visc;
}
}
}
i have tried to explain what each loop is doing incase you want any
clarification plz let me know if i can copy and paste it in your
compiler and could help me out in much more details it would be
grateful
Sachin Kanetkar