Hello Everyone
I am a beginner in this field and trying simple codes in LBM.I have wriiten the code but is not working and not giving suitable results. If some one can help me.
here is the code
#include<stdio.h>
#include"crun4.h"
void main()
{
FILE *a,ifp,fp3,fp4;
a=fopen(“lid.dat”,“w”);
ifp = fopen(“grid.dat”, “w”);
fp3=fopen(“uvel.dat”,“w”);
fp4 = fopen(“vvel.dat”, “w”);
fprintf(ifp,“ZONE I=200\t J=200”);
fprintf(ifp,"\n");
dx=1.0/nx;
dy=1.0/ny;
dt=dx;
tau=0.7;
nu=((2tau-1.0)dxdx)/(6dt);
u1=(nuRex);
x[0][0]=0.0;
y[0][0]=0.0;
x[nx][j]=1.0;
y[i][ny]=1.0;
/*----------------------------grid generation----------------------------*/
for(i=1;i<nx;i++)
{
for(j=1;j<ny;j++)
{
x[i][j]=i*dx;
y[i][j]=j*dy;
fprintf(ifp,"%lf\t %lf\n", x[i][j],y[i][j]);
//printf("grids are: %f\t%f\n",x[i],y[j]);
}
}
/*----------------------------initial condition----------------------------*/
for(i=0;i<=nx;i++)
{
for(j=0;j<ny;j++)
{
rho[i][j]=1.0;
ux[i][j]=0.0;
uy[i][j]=0.0;
for(k=0;k<9;k++)
{
f[i][j][k]=w_f[k]*rho[i][j];
fp[i][j][k]=w_f[k]*rho[i][j];
// //printf(“initialised PDF\n”);
}
}
}
for(i=1;i<nx;i++)
{
rho[i][ny]=1.0;
ux[i][ny]=u1;
uy[i][ny]=0.0;
// F2=0.01;
// for(k=0;k<9;k++)
// {
// F1=ex[k]0.1;
// feq=w_f[k]rho[i][ny](1.0+(3.0F1)+(4.5F1F1)-(1.5*F2));
// f[i][ny][k]=feq;
// //printf("%f\n",f[i][ny][k]);
// }
}
for(i=0;i<=nx;i++)
{
// rhon=f[i][ny][0]+f[i][ny][1]+f[i][ny][3]+2.0*((f[i][ny][2])+(f[i][ny][6])+(f[i][ny][5]));
F2=u1*u1;
for(k=0;k<9;k++)
{
F1=ex[k]*u1;
feq=w_f[k]*rho[i][ny]*(1.0+(3.0*F1)+(4.5*F1*F1)-(1.5*F2));
f[i][ny][k]=feq;
// //printf("%f\n",f[i][ny][k]);
}
/f[i][ny][4]=f[i][ny][2];
f[i][ny][8]=f[i][ny][6]+(rho[i][ny](ux[i][ny]/6.0));
f[i][ny][7]=f[i][ny][5]+(rho[i][ny](ux[i][ny]/6.0));/
}
/*-------------------------------main loop---------------------------------*/
error=10.0;
do
{
for(i=1;i<nx;i++)
{
for(j=1;j<ny;j++)
{
E2=ux[i][j]ux[i][j]+uy[i][j]uy[i][j];
for(k=0;k<9;k++)
{
E1=ex[k]ux[i][j]+ey[k]uy[i][j];
feq=w_f[k]rho[i][j](1.0+(3.0E1)+(4.5E1E1)-(1.5E2));
f[i][j][k]=fp[i][j][k]+((feq-fp[i][j][k])/tau);
//printf("%f\n",fp[i][j][k]);
}
}
}
/*-----------------------------Streaming step-----------------------------*/
for(j=0;j<=ny;j++)
{
for(i=nx;i>=1;i--)
{
fp[i][j][1]=f[i-1][j][1];
}
for(i=0;i<nx;i++)
{
fp[i][j][3]=f[i+1][j][3];
}
}
for(j=ny;j>=1;j--) /*At j=0,we have pdfs from Boundary conditions*/
{
for(i=0;i<=nx;i++)
{
fp[i][j][2]=f[i][j-1][2];
}
for(i=nx;i>0;i--)
{
fp[i][j][5]=f[i-1][j-1][5];
}
for(i=0;i<nx;i++)
{
fp[i][j][6]=f[i+1][j-1][6];
}
}
for(j=0;j<ny;j++) /*At j=ny,pdfs are defined*/
{
for(i=0;i<=nx;i++)
{
fp[i][j][4]=f[i][j+1][4];
}
for(i=0;i<nx;i++)
{
fp[i][j][7]=f[i+1][j+1][7];
}
for(i=nx;i>=1;i--)
{
fp[i][j][8]=f[i-1][j+1][8];
}
}
/*---------------------------------------Boundary Condition--------------------------------------*/
/*---------------------left wall boundary condition-----------------------*/
for(j=1;j<ny;j++)
{
f[0][j][1]=fp[0][j][3];
f[0][j][5]=fp[0][j][7];
f[0][j][8]=fp[0][j][6];
}
/*---------------------right wall boundary condition-----------------------*/
for(j=1;j<ny;j++)
{
f[nx][j][3]=fp[nx][j][1];
f[nx][j][6]=fp[nx][j][8];
f[nx][j][7]=fp[nx][j][5];
}
/*---------------------bottom wall boundary condition-----------------------*/
for(i=0;i<=nx;i++)
{
f[i][0][2]=fp[i][0][4];
f[i][0][5]=fp[i][0][7];
f[i][0][6]=fp[i][0][8];
}
/*-----------------------------Calculation of velocity and density-----------------------------*/
for(i=1;i<nx;i++)
{
for(j=1;j<ny;j++)
{
sum=0.0;
for(k=0;k<9;k++)
{
sum=sum+f[i][j][k];
}
rho[i][j]=sum;
// rho[i][j]=f[i][j][1]+f[i][j][2]+f[i][j][3]+f[i][j][4]+f[i][j][5]+f[i][j][6]+f[i][j][7]+f[i][j][8]+f[i][j][0];
// ux[i][j]=((f[i][j][1]+f[i][j][5]+f[i][j][8])-(f[i][j][3]+f[i][j][6]+f[i][j][7]))/rho[i][j];
// uy*/[i][j]=((f[i][j][2]+f[i][j][5]+f[i][j][6])-(f[i][j][4]+f[i][j][7]+f[i][j][8]))/rho[i][j];
}
}
for(i=1;i<nx;i++)
{
for(j=1;j<ny;j++)
{
uxsum=0.0;
uysum=0.0;
for(k=0;k<9;k++)
{
uxsum=uxsum+(f[i][j][k]*ex[k]);
uysum=uysum+(f[i][j][k]*ey[k]);
}
ux[i][j]=uxsum/rho[i][j];
uy[i][j]=uysum/rho[i][j];
}
}
/*-------------------------------------error checking condition------------------------------------------*/
e1=e2=0.0;
for(i=1;i<nx;i++)
{
for(j=1;j<ny;j++)
{
e1=e1+(sqrt((ux[i][j]-ux0[i][j])*(ux[i][j]-ux0[i][j])+(uy[i][j]-uy0[i][j])*(uy[i][j]-uy0[i][j])));
e2=e2+(sqrt((ux[i][j]*ux[i][j])+(uy[i][j]*uy[i][j])));
ux0[i][j]=ux[i][j];
uy0[i][j]=uy[i][j];
}
}
error=e1/e2;
if(timestep%1000==0)
{
printf("%d\t%0.13f\n",timestep,error);
}
timestep++;
// if(timestep>=100000) error=0.0;
}
while(error>=1e-8);
printf("%d\t%0.13f\n",timestep,error);
/---------------------------------------------output------------------------------------------------------/
fprintf(a,“Zone I=%d\tJ=%d \n”,nx+1,ny+1 );
for(j=0;j<=ny;j++)
{
for(i=0;i<=nx;i++)
{
fprintf(a,"%f\t%f\t%f\t %f\n",x[i],y[j],ux[i][j],uy[i][j]);
}
}
for(j=1;j<=ny;j++)
{
fprintf(fp3,"\n%f\t%f",y[100][j],ux[100][j]);
}
for(i=1;i<=nx;i++)
{
fprintf(fp4,"\n%f\t%f",x[i][100],uy[i][100]);
}
fclose(ifp);
fclose(a);
fclose(fp3);
fclose(fp4);
//getch();
}