[a_gogoi@yahoo.com]

i am trying a lid cavity problem using stream function-vorticity formulation. However the c code developed stops after 8 to 10 time steps . I am attaching the code and would be grateful for suggestions.
i have tried to run the program with different values of Array siz , time step but the program goes out of bounds.

# include <stdio.h>
# include <math.h>
# include <string.h>

# define SIZE 50

// assume id=jd ... no of divisions in x axis and y axis equal
main()
{

int i,j,k,l,m,n,id,maxiter; // re :reynolds no,id= grid size .. assume uniform mesh
double x[SIZE],y[SIZE],shi[SIZE][SIZE],zeta[SIZE][SIZE],u[SIZE][SIZE],v[SIZE][SIZE],a,b,c;
double maxstep,shires,re,h,t,dt,w;
//char chk;

FILE *fpshi,*fpzeta,*fpuv;

void zetabc( double shi[SIZE][SIZE],double zeta[SIZE][SIZE],int id,double h); // function decleration for zeta boundary conditions

//initialising the values of shi,zeta etc..
id=SIZE;
maxstep=1000;
maxiter=2000;
m=0;
;
w=1.8;// relaxation factor

re=500.00;
h=1.0/(id-1);
t=0.0;
dt=1.00/maxstep;

for(i=0;i<id;i++)
for(j=0;j<id;j++)
{
shi[i][j]=0.0;
zeta[i][j]=0.0;
u[i][j]=0.0;v[i][j]=0.0;
}
//setting boundary conditions for zeta
zetabc(shi,zeta,id,h);

for(i=0;i<id;i++)
{
x[i]=i*h;
y[i]=i*h;
}

//Starting Iteration process

for(k=0;k<maxstep;k++)
{
for(l=0;l<maxiter;l++)
{
// solving for stream function for internal mesh points
for(i=1;i<(id-1);i++)
{
for(j=1;j<(id-1);j++)
{
// residual of shi
shires=((0.25*(shi[i+1][j]+shi[i-1][j]+shi[i][j-1]+shi[i][j+1]+((h*h)*zeta[i][j])))-shi[i][j]);
printf("shires %d \t %d \t %4.15f \n",i,j,shires);
// SOR
shi[i][j]=shi[i][j]+(w*shires);
}
}




}

// set boundary conditions for vorticity
zetabc(shi,zeta,id,h);


// computing vorcitity equation for internal mesh points
for(i=1;i<(id-1);i++)
{ for(j=1;j<(id-1);j++)
{
a=(1/(re*h*h))*(zeta[i+1][j]+zeta[i-1][j]+zeta[i][j+1]+zeta[i][j-1]-(4.0*zeta[i][j]));
b=(0.25/(h*h))*((shi[i+1][j]-shi[i-1][j])*((zeta[i][j+1]-zeta[i][j-1])));
c=(0.25/(h*h))*((shi[i][j+1]-shi[i][j-1])*(zeta[i+1][j]-zeta[i-1][j]));

zeta[i][j]=zeta[i][j]+ t*(a+b-c);
printf("zeta is %f\n",zeta[i][j]);

}
}

printf("out of loop\n");
t=t+dt; // increment time step

//computing data for u and v form stream function
for(i=1;i<(id-1);i++)
{ for(j=1;j<(id-1);j++)
{

u[i][j]=(shi[i][j+1]-shi[i][j-1])/(2 * h);
v[i][j]=-1.0*(shi[i+1][j]-shi[i-1][j])/(2 * h);
u[i][id-1]=1.00;

}

}
printf("after uv computation\n");
// writing values of shi,

fpshi=fopen("shi.txt","w");
fpzeta=fopen("zeta.txt","w");
fpuv=fopen("uv.txt","w");
printf("After opening files\n");

for(i=0;i<id;i++)
for(j=0;j<id;j++)
{
fprintf(fpshi, "%2.5f \t %2.5f \t %20.10f \n",x[i],y[j],shi[i][j]);
fprintf(fpzeta,"%2.5f \t %2.5f \t %20.10f \n",x[i],y[j],zeta[i][j]);
fprintf(fpuv, "%2.5f \t %2.5f \t %20.10f \t %20.10f\n",x[i],y[j],u[i][j],v[i][j]);
printf("write data %d %d\n",i,j);
}

fclose(fpshi);
fclose(fpzeta);
fclose(fpuv);

}
}

// setting boundary conditions for vorticity zeta
void zetabc( double shi[SIZE][SIZE], double zeta[SIZE][SIZE],int id,double h)
{
int i,j;
printf("in zeta BC\n");

//top and bottom wall
for(i=0;i<=(id-1);i++)
{

zeta[i][0]=(2/(h*h))*(-shi[i][1]);
zeta[i][id-1]=((2/(h * h))*(-shi[i][id-2]))-(2.0/h);
printf("zeta top and bottom \t %2.10f \t %2.10f\n",zeta[i][0],zeta[i][id-1]);

}
printf("finish top n bottom wall \n");

// left and right wall
for (i=1;i<=(id-2);i++)
{
zeta[0][i]=(2/(h*h))*(-shi[1][i]);
zeta[id-1][i]=(2/(h*h))*(-shi[id-2][i]);
printf("zeta left and right \t %2.10f \t %2.10f\n",zeta[0][i],zeta[id-1][i]);
}
return;
}

[Alisha]

hi Arunav,
I guess you have not set the initial condition of vorticity properly.

Though I don't code with C ( was struggling to understand it ), it seems like you have set vorticity 0 everywhere, which is not the case... in fortran I would write it like this :

DO i=0,n
DO j=0,n

IF (j==1) THEN
vorticity(i,j) = u_lid/dy
ELSE
vorticity(i,j) = 0
END IF

END DO
END DO

which means at moving wall the vorticity wont be 0
Tell me if it solves your problem...
and first try with no relaxation i.e. w = 1, it will work too... take Re = 100 for the first attempt...

Additionally I am sending the link of a pdf which I found helpful during writing stream fn vorticty code... Happy coding

[a_gogoi@yahoo.com]

Dear Alisha,
i start my expressing my gratitude for your reply.
I have done the following in my code.
Before starting ,in C unlike FORTRAN for an array of a[id], the arrays start with a[0] and end with a[id-1].
Also in my code i have used zeta for vorticity and shi for stream function.
void zetabc(...) is function for setting the boundary conditions for zeta.
hence the problem of zeta is i think taken care of .
I have also tried your suggestions my i am unable to get the results.

[Alisha]

I have a small fortran code, which you can use if u wish. Just give me your mail id

[a_gogoi@yahoo.com]

my email id is a_gogoi@yahoo.com