I write a c++ program to model free convection heat transfer using LBM ( of course!).
In order to add the force term to the collision stage, I noticed that there are several ways to do so.
One of them is by adding to f an extra df(i) of the form: 3w(i)cF/cs^2.
I tried to develop this expression and couldn’t.

Here I present my expression for D2Q9:

df(1)c(1)+…+ df(9)c(9)=rowu (row is the density and u is the macroscopic velocity).
rowu=rowadt=F*t=F (assuming dt equal 1 at LB units).

df(i)c(i)=w(i)*F
df(i)=w(i)*dot(F,c(I))/c(i)^2

c(i) equal 1 at LB units, so:

dfi(i)= w(i)*dot(F,c(I))

I will really appreciate any comment regarding the way I developed I would which
you could tell me please.

My apologies; I read the ‘c’ in your first expression as c_i. This is the basic form of a mass conserving acceleration term. The force as you have it (R_i~W_i*constant) will change the mass but not the momentum.