Hello everyone. I am an undergraduate student, learning lbm for my btech project. I started with sukop and succi books. I tried poiseulle flow under gravity, but i cant understand where the fault is. Please can any one help me in this or any other further reference material would be helpful. Here is my code in c
#include<stdio.h>
main()
{
double Re, umax, uavg, g, v, tou,temp,z;
int LX, LY,a,i,j,t,ip,in,jp,jn,x,l;
double d = 25;
l=50;
double u[50];
Re = 10.0;
umax = 0.1;
uavg = (2.0/3.0)umax;
LX = 51;
LY = 600;
v = (uavgLX)/Re;
tou = ((6.0v)+1.0)/2.0;
printf("%lf",tou);
g = 3.0uavgv(2.0/LX)*(2.0/LX);
printf("%lf" , g);
double u_x[LX][LY];
double u_y[LX][LY];
double ueq_x[LX][LY];
double ueq_y[LX][LY];
double rho[LX][LY];
double feq[LX][LY][9];
double f[LX][LY][9];
double ftemp[LX][LY][9];
int solid_node[LX][LY];
double uxsq[LX][LY];
double uysq[LX][LY];
double usq[LX][LY];
double w1 = 4.0/9.0;
double w2 = 1.0/9.0;
double w3 = 1.0/36.0;
double ex[9];
ex[0] = 0.0;
ex[1] = 1.0;
ex[2] = 0.0;
ex[3] = -1.0;
ex[4] = 0.0;
ex[5] = 1.0;
ex[6] = -1.0;
ex[7] = -1.0;
ex[8] = 1.0;
double ey[9];
ey[0] = 0.0;
ey[1] = 0.0;
ey[2] = 1.0;
ey[3] = 0.0;
ey[4] = -1.0;
ey[5] = 1.0;
ey[6] = 1.0;
ey[7] = -1.0;
ey[8] = -1.0;
for(j=0; j<LY; j++)
{
for(i=0;i<LX;i++)
{
if(i==0 || i==(LX-1))
{
solid_node[i][j] = 1;
}
else
{
solid_node[i][j] = 0;
}
}
}
//initialization for rho u and f
for(i=0; i < LX; i++)
{
for(j=0; j < LY; j++)
{
rho[i][j] = 1.0;
u_x[i][j] = 0.0;
u_y[i][j] = 0.0;
f[i][j][0] = w1;
f[i][j][1] = w2;
f[i][j][3] = w2;
f[i][j][4] = w2;
f[i][j][2] = w2;
f[i][j][5] = w3;
f[i][j][6] = w3;
f[i][j][7] = w3;
f[i][j][8] = w3;
}
}
// time loop
t = 0;
while(t<30000)
{
t++;
//macroscopic velocities
for(i=0; i < LX; i++)
{
for(j=0; j < LY; j++)
{
rho[i][j] = 0.0;
u_x[i][j] = 0.0;
u_y[i][j] = 0.0;
if(!solid_node[i][j])
{
for(a=0;a<9;a++)
{
rho[i][j] = rho[i][j] + f[i][j][a];
u_x[i][j] = u_x[i][j] + ex[a]*f[i][j][a];
u_y[i][j] = u_y[i][j] + ey[a]*f[i][j][a];
}
u_x[i][j] = u_x[i][j]/rho[i][j];
u_y[i][j] = u_y[i][j]/rho[i][j];
}
}
}
//equilibrium distribution functions
for(i=0; i < LX; i++)
{
for(j=0; j < LY; j++)
{
if(!solid_node[i][j])
{
ueq_x[i][j] = u_x[i][j];
ueq_y[i][j] = u_y[i][j] + tou*g;
uxsq[i][j] = ueq_x[i][j]*ueq_x[i][j];
uysq[i][j] = ueq_y[i][j]*ueq_y[i][j];
feq[i][j][0] = w1*rho[i][j]*(1 - (3.0/2.0)*(uxsq[i][j]+uysq[i][j]));
feq[i][j][1] = w2*rho[i][j]*(1+3.0*ueq_x[i][j]+(9.0/2.0)*ueq_x[i][j]*ueq_x[i][j]-(3.0/2.0)*(uxsq[i][j]+uysq[i][j]));
feq[i][j][2] = w2*rho[i][j]*(1+3.0*ueq_y[i][j]+(9.0/2.0)*ueq_y[i][j]*ueq_y[i][j]-(3.0/2.0)*(uxsq[i][j]+uysq[i][j]));
feq[i][j][3] = w2*rho[i][j]*(1-3.0*ueq_x[i][j]+(9.0/2.0)*ueq_x[i][j]*ueq_x[i][j]-(3.0/2.0)*(uxsq[i][j]+uysq[i][j]));
feq[i][j][4] = w2*rho[i][j]*(1-3.0*ueq_y[i][j]+(9.0/2.0)*ueq_y[i][j]*ueq_y[i][j]-(3.0/2.0)*(uxsq[i][j]+uysq[i][j]));
feq[i][j][5] = w3*rho[i][j]*(1+3.0*(ueq_x[i][j]+ueq_y[i][j])+(9.0/2.0)*(ueq_x[i][j]+ueq_y[i][j])*(ueq_x[i][j]+ueq_y[i][j])-(3.0/2.0)*(uxsq[i][j]+uysq[i][j]));
feq[i][j][6] = w3*rho[i][j]*(1+3.0*(-ueq_x[i][j]+ueq_y[i][j])+(9.0/2.0)*(-ueq_x[i][j]+ueq_y[i][j])*(-ueq_x[i][j]+ueq_y[i][j])-(3.0/2.0)*(uxsq[i][j]+uysq[i][j]));
feq[i][j][7] = w3*rho[i][j]*(1+3.0*(-ueq_x[i][j]-ueq_y[i][j])+(9.0/2.0)*(-ueq_x[i][j]-ueq_y[i][j])*(-ueq_x[i][j]-ueq_y[i][j])-(3.0/2.0)*(uxsq[i][j]+uysq[i][j]));
feq[i][j][8] = w3*rho[i][j]*(1+3.0*(ueq_x[i][j]-ueq_y[i][j])+(9.0/2.0)*(ueq_x[i][j]-ueq_y[i][j])*(ueq_x[i][j]-ueq_y[i][j])-(3.0/2.0)*(uxsq[i][j]+uysq[i][j]));
}
}
}
// streaming
for(i = 0; i < LX; i++)
{
for(j=0; j < LY; j++)
{
ip = ( i < LX-1)?(i+1):(0);
in = ( i > 0)?(i-1):(LX-1);
jp = ( j < LY-1)?(j+1):(0);
jn = ( j > 0)?(j-1):(LY-1);
ftemp[i][j][0] = f[i][j][0];
ftemp[ip][j][1] = f[i][j][1];
ftemp[i][jp][2] = f[i][j][2];
ftemp[in][j][3] = f[i][j][3];
ftemp[i][jn][4] = f[i][j][4];
ftemp[ip][jp][5] = f[i][j][5];
ftemp[in][jp][6] = f[i][j][6];
ftemp[in][jn][7] = f[i][j][7];
ftemp[ip][jn][8] = f[i][j][8];
}
}
//collision
for(i = 0; i< LX; i++)
{
for(j=0; j < LY; j++)
{
if(!solid_node[i][j])
{
for(a=0;a<9;a++)
{
f[i][j][a] = ftemp[i][j][a]-((ftemp[i][j][a]-feq[i][j][a])/tou);
}
}
else
{
temp = f[i][j][1];
f[i][j][1] = f[i][j][3];
f[i][j][3] = temp;
temp = f[i][j][2];
f[i][j][2] = f[i][j][4];
f[i][j][4] = temp;
temp = f[i][j][5];
f[i][j][5] = f[i][j][7];
f[i][j][7] = temp;
temp = f[i][j][6];
f[i][j][6] = f[i][j][8];
f[i][j][8] = temp;
}
}
}
}
//output
for(i=0; i < LX; i++)
{
for(j=0; j < LY; j++)
{
if(!solid_node[i][j])
{
rho[i][j] = 0.0;
u_x[i][j] = 0.0;
u_y[i][j] = 0.0;
for(a=0;a<9;a++)
{
rho[i][j] = rho[i][j] + f[i][j][a];
u_x[i][j] = u_x[i][j] + ex[a]*f[i][j][a];
u_y[i][j] = u_y[i][j] + ey[a]*f[i][j][a];
}
u_x[i][j] = u_x[i][j]/rho[i][j];
u_y[i][j] = u_y[i][j]/rho[i][j];
printf("(%d,%d).....(%lf,%lf)\n ",i,j,u_x[i][j],u_y[i][j]);
}
}
}
}