Regularized BC and half-way bounce-back BC for corners


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: (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];

Hello Bart,

Why do you say Regularized (REG) method isn’t valid for corners? I believe the REG is a local method hence you should be able to use it at corners. (perhaps some special reasoning must be used for the populations belonging to buried links as in the Zou/He approach).

As the external forcing term in LB enters Chapman-Enskog expansion at first order in \epsilon I would say that it is affecting the f(neq) term. Yet, as your force models gravity (i.e. is constant both in space and time) I think (though I am not 100% sure) that the error you introduce, by not considering the force in f(neq) formulation, will end up vanishing at macroscopic level.



Thank you for your answer. I have done some experimenting and read Zou/He’s paper among others.

Let’s take the case of a bottom-left internal corner:

6   2   5
  \ | / 
3 -- -- 1
  / | \
7   4   8

f3, f4, and f7 stream from the fluid and are valid, the rest stream from the wall and are therefore unknown. Now, it is already not straightforward to determine the density (assuming velocity is known). But in [1] a formula is given as follows:
rho = 36 * f7
Now it is possible to compute the equilibria of all the f’s. RLBM assumes we can bounce back off-equilibria, so
neq1 = neq3
neq2 = neq4
neq5 = neq7
With this info we can compute (a non-regularized) f1, f2, and f5 by summing their corresponding eq and neq.

Then we still need to determine f6 and f8. In Zou/He’s paper, they compute f6 and f8 as follows:
f6 = f8 = 0.5 * (rho - (f1 + f2 + f3 + f4 + f5 + f7))

Then we compute neq6 and neq8, so all neq’s are known. Then we execute the regularization procedure by compute Pi^neq, etc. and determine a regularized version of all the f’s.

Doing this, the artifacts near the corners are less, but they are still present. Also, the second issue that the total density of the fluid keeps decreasing is still present. Is RLBM BC by the way mass-conserving are not?

If you (or anybody) else have suggestion, please let me know.

The gravity force term indeed affects the neq-part. But this gravity force term is applied to the f’s in the collision-operator, and the f’s stream to neighbors and therefore the f’s with gravity force will still be present in handling the boundaries. So I would think it is odd that we would include the force term again in handling the boundary conditions. We are also not including the gravity force term again if we perform the half-way no-slip BC (nor in full-way no-slip BC).

But having said all this. Why wouldn’t it work to have the combination of half-way bounce-back BC for corners and regularized BC for straight boundaries? If have also seen people use Skordos finite-difference BC for corners and regularized BC for straight boudnaries.

[1] Joris C. G. Verschaeve. Analysis of the lattice Boltzmann Bhatnagar-Gross-Krook no-slip boundary condition: Ways to improve accuracy and stability. Phys. Rev. E 80, 036703 (2009)

Hello Bart,

All questions you have raised are quite interesting indeed. Yet I do not know if I can give you the correct answers for all of them. Here are my attempts:

  1. Regarding the corner issue.

All steps you use in your procedure seem similar to those I employ. All, except the rho computation. I extrapolate the rho corner value from inside. The rho formula you use belongs to a mass-conserving scheme. Zou/He is not mass conserving (and I suspect REG isn’t neither). For instance, in Verschaeve’s mass-conserving scheme the non-equilibrium populations are not bounced-back but fneq(i) = -fneq(iopp). Hence you use a mix of schemes.Yet, I do not know how harmful that can be. Furthermore, and if you are simulating incompressible flow I wouldn’t worry on the mass conservation issue. In the incompressible regime mass stands for pressure and what counts here are pressure differences. The pressure (mass) value by itself is meaningless.

  1. Regarding the body force.

It is true that F (the force) affects f (the distributions) in the collision step. However, in order to account for the discrete lattice effects, or in other words, in order to recover the correct forcing Navier-Stokes equations, with the same second order accuracy of the unforced LB method, the force must be considered in the computation of the velocity. Therefore, the velocity is given not only by the f’s first order velocity moment but it is also affected by F.
Note for instance Zou/He model with a force per unit mass g used to model a resting wall. (I use the incompressible model! Therefore, rho0 is constant in j=rho0*u)

    for x=1:LX



For a corner:

    % Bottom left (inlet)
    x=1; y=1;


I hope you can have an idea from these examples.

Good luck.



I did some more experimenting with the corner:
Extrapolating the density from one or more neighbors gives me more or less the same artifacts near the corners as when I used halfway bounce-back for the corners. So this is worse than using the formula in [1]. The stability of the simulation is also worse than using halfway bounce-back for all boundaries.

Getting rid of regularized BC for corners and using the approach of [1] for corners, but still use regularized BC for straight walls, removes the artifacts near the corners. The stability is more or less equal as using halfway bounce-back for all boundaries. The decreasing rho is still present though, but as you said it, this might be inherent of the regularized BC. In [1] it is reasoned that neq = 0 for the f’s at an internal corner, so the f’s receive the equilibrium distribution functions (I did not check all his formulas since I no longer have the time to sort these things out in detail)

As for the body force. There are multiple formulas for implementing a gravity force in the collision operator, for some of these formulas also need you to correct (some of the) moments accordingly. But not all formulas have this requirement, e.g. Guo’s formula does not require this, and as far as I know the formula I currently use, Fi = (ei . gravity) * rho * wi / cs2, also does not require this.