Hello,
I am currently working on regularized LBM and now regularized LBM with half-way bounce-back BC is finished and working properly, I also wanted to have regularized boundary conditions. However, here a problem appears in combination with a gravity force.
The situation we have is a 2D box with the D2Q9-model and the force term Fi = (ei . gravity) * rho * wi / cs2. Since regularized BC is not appropriate for corner nodes, since necessary information is missing, I decided to use half-way bounce back for the corners. So all all fluid nodes with a straight boundary neighbor are handled with regularized BC and all fluid nodes with with a corner boundary neighbor are handled with half-way bounce back BC. I obtain some weird results, seemingly coming from the corners. The regularized BC keeps the velocity zero, but near the corners the velocity is non-zero (I made a screenshot of the situation: http://home.kpn.nl/riji99ou/screenshot.png (black nodes are boundary)). The second issue is that the total density of the fluid keeps decreasing.
If I run the lid-driven cavity with the same conditions as above (but no gravity of course), there are no apparent errors. So the questions pops up: do I have to include the gravity force in regularized BC somewhere, or handle the gravity force somehow? (would seem strange though since I already have a gravity force included in the collision operator)
I compared my code to the C-version in the section with sample code, this version also includes regularized BC, and I could not find any difference in output of rho and Pi^neq. Although I have confidence in the code, I’ll post the code for the upper boundary below (perhaps I am missing something).
Any help or hints would greatly be appreciated.
//Velocity vectors
// 0 1 2 3 4 5 6 7 8
int ex[9] = {0, 1, 0, -1, 0, 1, -1, -1, 1};
int ey[9] = {0, 0, 1, 0, -1, 1, 1, -1, -1};
void LBMsim::computeUpperRegularized(int x, int y) {
//Compute rho
double sumUpper = f[next][x][y][2] + f[next][x][y][5] + f[next][x][y][6];
double sumMiddle = f[next][x][y][0] + f[next][x][y][1] + f[next][x][y][3];
double uyCell = uy[x][y];
rho[x][y] = 1.0 / (1.0 + uyCell) * (2.0 * sumUpper + sumMiddle);
double eq[9];
double neq[9];
double piNeq_xx = 0.0, piNeq_yy = 0.0, piNeq_xy = 0.0;
for (int i = 0; i < 9; ++i) {
double fi = f[next][x][y][i];
//Compute equilibrium distribution function
eq[i] = fEq(i, rho[x][y], ux[x][y], uy[x][y]);
//Compute off-equilibrium distribution function
neq[i] = fi - eq[i];
}
//Bounce-back off-equilibrium parts
neq[4] = neq[invDir[4]];
neq[7] = neq[invDir[7]];
neq[8] = neq[invDir[8]];
for (int i = 0; i < 9; ++i) {
//Compute off-equilibrium momentum flux tensor
piNeq_xx += neq[i] * ex[i] * ex[i];
piNeq_yy += neq[i] * ey[i] * ey[i];
piNeq_xy += neq[i] * ex[i] * ey[i];
}
double f1[9];
for (int i = 0; i < 9; ++i) {
//Compute first-order term of distribution function
double Qixx = ex[i] * ex[i] - cs2;
double Qiyy = ey[i] * ey[i] - cs2;
double Qixy = ex[i] * ey[i];
f1[i] = w[i] / (2.0 * cs2 * cs2) * (Qixx * piNeq_xx + Qiyy * piNeq_yy + 2.0 * Qixy * piNeq_xy);
//Replace f by regularized counterpart
f[next][x][y][i] = eq[i] + f1[i];
}
}