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.102*10^-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.5omega)*ω (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.5F(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.5F(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