# Gravity Force in Lattice Boltzmann

Dear All,

I have problem with implementation of body forces in the LBM. Considering the Poiseuille flow in a slit driven by gravity, I don’t know how to incorporate gravitational acceleration into the Model. Based on the Sukop, u^eq , should be used in the computation of equilibrium distribution function:
u^eq=u+∆u=u+ (τ(ρg))/ρ
The gravitational force have component only in the vertical direction, so whether the above equation have to be applied only to the velocity along the slit or not?
Does the algorithm of macro properties calculation change or not?

Any help is appreciated,
Regards…

Hi Narges,

I think that the articles by Guo (2002) and Kupershtokh (2004) are indeed good references to understand about forcing in LBM. Please post any questions once you are through with these articles.

Best,
Timm

Dear Timm,

I read the paper by Guo, and chose the first method of body force implementation. The problem to be solved is (as reported by Sukop):
yDim = 11 lu
Re = 4.4
Thaw = 1/omega = 1.0
Kinematic viscosity = 0.166667
Maximum velocity = 0.1
Gravitational acceleration = g = 1.10210^-3
the collision term is changed to :
if (image(y,x)≠ wall)
f(y,x,i) = (1.0d0 - omega) * f(y,x,i) + omega * feq(y,x,i) + Force(i)
where,
Force(i) = (1- 0.5
omega)ω (i){[ 3*(e(i,0)-u(y,x,0)) + 9*e(i,0)*u(y,x,0) + e(i,1)*u(y,x,1))e(i,0) )]F(0) + [ 3(e(i,1)-u(y,x,1)) + 9e(i,0)*u(y,x,0) + e(i,1)*u(y,x,1))e(i,1) )]F(1)}
In this case, F(0)= ρ(y,x)g and F(1)=0.
The other change that must be applied is the macro properties calculation as follows:
u(y,x,0) =[f(y,x,1) - f(y,x,3) + f(y,x,5) - f(y,x,6) - f(y,x,7) + f(y,x,8) + 0.5
F(0)] / rho(y,x)
u(y,x,1) =[f(y,x,2) - f(y,x,4) + f(y,x,5) + f(y,x,6) - f(y,x,7) - f(y,x,8) + 0.5
F(1)] / rho(y,x)
Would you please tell me whether the above formulas are correct or not?
And the algorithm shown below is correct or not?
do tStep = 1,tMax
CALL wall_boundaries(f,image)
CALL compute_Macros(f,rho,u,uSqr)
CALL compute_Feq(fEq,rho,u,uSqr)
CALL collide(f,fEq,omega,image,rho,u)
CALL stream(f) ! with regard to periodic boundary conditions
new=sum(u(:,:,0))
t=(new-old)/old
print
,tstep,abs(t)
if(abs(t).lt.eps) then
exit
end if
old=new
end do
I highly appreciate your considerations before,
Sincerely,
Narges

Hi Narges,

as you mentioned from Sukop book, the gravitational acceleration is incorporated in a velocity term and it would be effect on the fluid only in the gravity direction , so you should apply the correct velocity (u+∆u) to obtain [f][/eq]. as well as, use bounce back BC along the state walls (left and right walls) and periodic boundary in the flow (top and bottom boundary). to get simple code which includes external force go to Mike Sukop website.

Best Regards,