Hi Dear Jonny,
Dear friends,
I am a begginer to the LBm. i want to implement Zuo/He bounce back scheme to my chamber problem. my problem description is given below
inlet outlet
__P1 P2_______P3 P4 ___
| |
| |
| wall | wall
| |
|______________________|
wall
my doubts are
- how to implement outflow boundary condition?
2.shall i have to implement the “concept of corner nodes” near the edge of inlet (2 points P1 and P2) and outlet ( 2points P2 and P4)?
I have written my code in C++. the boundary conditions are as follows. my distribution function notation is from 0 to 8
void setBoundary()
{
// ... left boundary condition
// ... Zou/He scheme
for( int j = 1; j < (NY-1); ++j )
{
Fout[0][j][1] = Fout[1][j][3];
Fout[0][j][5] = Fout[1][j+1][7];
Fout[0][j][8] = Fout[1][j-1][6];
//Left Corner down node
Fout[0][0][1] = Fout[1][1][3];
Fout[0][0][5] = Fout[1][1][7];
Fout[0][0][2] = Fout[1][1][4];
Fout[0][0][6] = 0.5*(rho[1][1]-(Fout[1][1][0]+Fout[1][1][1]+Fout[1][1][2]+Fout[1][1][3]+Fout[1][1][4]+Fout[1][1][5]+Fout[1][1][7]));
Fout[0][0][8] = Fout[0][0][6];
//left corner up side node
//rho[1][NY-2]= Fout[1][NY-2][0]+Fout[1][NY-2][2]+Fout[1][NY-2][3]+Fout[1][NY-2][6]+Fout[1][NY-2][7]+Fout[1][NY-2][5]+Fout[1][(NY-2)-1][2]+Fout[(1+1)][NY-2][3]+Fout[(1+1)][(NY-2)-1][6];
Fout[0][NY-1][1] = Fout[1][NY-2][3];
Fout[0][NY-1][8] = Fout[1][NY-2][6];
Fout[0][NY-1][4] = Fout[1][NY-2][2];
Fout[0][NY-1][5] = 0.5*(rho[1][1]-(Fout[1][NY-2][0]+Fout[1][NY-2][1]+Fout[1][NY-2][2]+Fout[1][NY-2][3]+Fout[1][NY-2][4]+Fout[1][NY-2][6]+Fout[1][NY-2][8]));
Fout[0][NY-1][7] = Fout[0][NY-1][5];
/*rho[j][0] = rho[j][1];
u_x[j][0] = uw;
u_y[j][0] = uw;*/
// right bounday as wall
Fout[NX-1][j][3] = Fout[NX-2][j][1];
Fout[NX-1][j][7] = Fout[NX-2][j-1][5];
Fout[NX-1][j][6] = Fout[NX-2][j+1][8];
//right corner down node
Fout[NX-1][0][6] = Fout[NX-2][1][8];
Fout[NX-1][0][2] = Fout[NX-2][1][4];
Fout[NX-1][0][3] = Fout[NX-2][1][1];
Fout[NX-1][0][5] = 0.5*(rho[NX-2][1]-(Fout[NX-2][1][0]+Fout[NX-2][1][1]+Fout[NX-2][1][2]+Fout[NX-2][1][3]+Fout[NX-2][1][4]+Fout[NX-2][1][6]+Fout[NX-2][1][8]));
Fout[NX-1][0][7] = Fout[NX-1][0][5];
//right corner up node
Fout[NX-1][NY-1][3] = Fout[NX-2][NY-2][1];
Fout[NX-1][NY-1][7] = Fout[NX-2][NY-2][5];
Fout[NX-1][NY-1][4] = Fout[NX-2][NY-2][2];
Fout[NX-1][NY-1][6] = 0.5*(rho[NX-2][NY-2]-(Fout[NX-2][NY-2][0]+Fout[NX-2][NY-2][1]+Fout[NX-2][NY-2][2]+Fout[NX-2][NY-2][3]+Fout[NX-2][NY-2][4]+Fout[NX-2][NY-2][5]+Fout[NX-2][NY-2][7]));
Fout[NX-1][NY-1][8] = Fout[NX-1][NY-1][6];
}
}
for(int i = 1; i < (NX-1); ++i)
{
//bottom bounday condition
Fout[i][0][2] = Fout[i][1][4];
Fout[i][0][5] = Fout[i+1][1][7];
Fout[i][0][6] = Fout[i-1][1][8];
// top bounday definition
if (i > 0 & (i<= NX/5))
{
Fout[i][NY-1][4] = Fout[i][NY-2][2];
Fout[i][NY-1][7] = Fout[i-1][NY-2][5];
Fout[i][NY-1][8] = Fout[i+1][NY-2][6];
}
else if ((i > (NX/5)) & (i <= (0.3*NX)))
{
u_x[i][NY-1] = 0.0;
u_y[i][NY-1] = -U;
rho[i][NY-1] = (1.0/(1-u_y[i][NY-1]))*((f[i][NY-1][0]+f[i][NY-1][1]+f[i][NY-1][3]) + 2*(f[i][NY-1][2]+f[i][NY-1][6]+f[i][NY-1][5]));
Fout[i][NY-1][4] = Fout[i][NY-1][2] + (2/3)*(rho[i][NY-1]*u_y[i][NY-1]);
Fout[i][NY-1][7] = Fout[i][NY-1][5] - 0.5*(Fout[i][NY-1][3]-Fout[i][NY-1][1]) + 0.5*(rho[i][NY-1]*u_x[i][NY-1]) + (1/6)*(rho[i][NY-1]*u_y[i][NY-1]);
Fout[i][NY-1][8] = Fout[i][NY-1][6] + 0.5*(Fout[i][NY-1][3]-Fout[i][NY-1][1]) - 0.5*(rho[i][NY-1]*u_x[i][NY-1]) + (1/6)*(rho[i][NY-1]*u_y[i][NY-1]);
}
else if ((i > (0.3*NX)) & (i <= (0.7*NX)))
{
Fout[i][NY-1][4] = Fout[i][NY-2][2];
Fout[i][NY-1][7] = Fout[i-1][NY-2][5];
Fout[i][NY-1][8] = Fout[i+1][NY-2][6];
}
else if((i > (0.7*NX)) & (i <= (4*NX/5)))
{
rho[i][NY-1] = rho[i][NY-2]; [NY-2][6]+Fout[i-1][NY-2][1]+Fout[i+1][NY-2][3]
u_x[i][NY-1] = -1.0 + (1.0/rho[i][NY-1])*((f[i][NY-1][0]+f[i][NY-1][1]+f[i][NY-1][3]) + 2*(f[i][NY-1][2]+f[i][NY-1][6]+f[i][NY-1][5])); u_y[i][NY-1] = 0.0; //u_y[i][NY-2];
Fout[i][NY-1][2] = Fout[i][NY-1][4]+ (2/3)*(rho[i][NY-1]*u_y[i][NY-1]);
Fout[i][NY-1][5] = Fout[i][NY-1][7] - 0.5*(Fout[i][NY-1][3]-Fout[i][NY-1][1]) + 0.5*(rho[i][NY-1]*u_x[i][NY-1]) + (1/6)*(rho[i][NY-1]*u_y[i][NY-1]);
Fout[i][NY-1][6] = Fout[i][NY-1][8] + 0.5*(Fout[i][NY-1][3]-Fout[i][NY-1][1]) - 0.5*(rho[i][NY-1]*u_x[i][NY-1]) + (1/6)*(rho[i][NY-1]*u_y[i][NY-1]);
}
else
{
Fout[i][NY-1][4] = Fout[i][NY-2][2];
Fout[i][NY-1][7] = Fout[i-1][NY-2][5];
Fout[i][NY-1][8] = Fout[i+1][NY-2][6];
}
}
}
Could you please have a look at my code and tell where i went wrong?
my density is increasing to very large number as iteration number increases. Please help me.
Thanks
winter