In 3D15Q, I want to add the gravity term as a body force.
According to J.M.Buick and C.A. Greated, Gravity in a lattice Boltzmann model(PHYSICAL REVIEW vol61,5),
I have to add [D/(bc^2)Fe] to the distribution function.

F = [ rho(at each node) * g_lattice] or [ rho(at each direction of a node) * g_lattice] ??

It seems to me that none of the two formulas does an appropriate job. The inclusion of a body force in lattice Boltzmann has been discussed in-depth in a paper by Guo e.a. Through a Chapman-Enskog expansion, the paper shows that a simple force-term, linear to g_lattice, gives raise to unwanted error terms. To get rid of them, an additional force-velocity coupling must be implemented. The appropriate formula can be found in the paper.

You are right, I am being too categorical with my answer: thank you for pointing this out. If the force field is space and time independent, YY’s force term is sufficient. Indeed, the error terms found in Guo’s Chapman-Enskog analysis are multiples of the time- resp. space-derivatives of the force. And… hmmm… now that I check the original post, I see that it is about gravity, which tends to be space and time independent.

So, the right way to implement the force term in this case is as follows. After collision, add the value F_i to each particle population:

F_i = rho t_i / c_s^2 c_i*g,

where rho is the density, t_i the lattice weight (e.g. 4/9, 1/9, 1/36 for D2Q9), c_s^2 the squared speed of sound (e.g. 1/3 for all standard lattices, including D3Q15), c_i the unitary lattice velocity ( e.g. (1,1) ) and g the force.

Thank you for your reply.
In my case, gravity is constant. Certainly, according to Guo, Buick’s method is not accurate.
But it may be enough.
Anyway, you mean rho is density at node , not at each direction of the node.
Is it right ?

Sukop’s implementation is not correct. I would stick to the articles by Guo or Khupershtokh (New method of incorporating a body force term into the lattice Boltzmann equation, Proc. 5th International EHD Workshop 2004).